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Abs  trac  t 


Using  a  self  adjoint  form  of  the  transport  equation 

expressed  as  a  variational  integral,  finite  element  equations 

for  the  one  dimensional,  one  speed,  homogeneous , time 

independent  transport  equation  in  slab  geometry  were  derived 

0  1 

and  encoded  in  Fortran  77.  The  accuracy  of  G  and  C  continuous 

fits  was  compared  against  an  analytical  solution  for  the  case 

0 

of  noscatter.  It  was  found  that  the  C  fits  require  an 

1 

excessive  amount  of  mesh  refinement.  The  C  fit  is  very 
accurate,  and  does  not  appear  to  be  computationally  excessive. 

The  finite  element  results  were  then  compared,  for  the  case 
of  isotropic  scatter,  to  a  legendre  polynomial  solution,  and 
the  results  of  a  recently  developed  code  known  as  Ln .  The 
methods  accuracy  was  sufficiently  verified  with  Inexact 
scattering  term  evaluation.  A  technique  of  exact  scattering 
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The  purpose  of  this  study  was  to  continue  the  work  of  a 

previous  graduate  student,  (A.D.  Goff  GNE  84M)  and  demonstrate 

that  a  finite  element  solution  of  the  transport  equation  would 

work.  Using  a  self  adjoint  form  of  the  transport  equation 

expressed  as  a  variational  integral,  finite  element  equations 
0  1 

with  C  and  C  continuity  were  derived,  encoded  and  compared  to 
a  spherical  harmonic  solution  over  a  test  case  domain. 
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being  used. 

-  Vector  of  finite  element  interpolating  nodes 

-  Vector  of  product  &  (p 

-  Vector  of  product  (p  Cp 


(  1 , 0 , 0  )  ,  (p  -  ^  a  t  aode  2,  and  ^  at  node  3.  This  fit  has  C 
continuity  since  flux  is  continuous  along  element  interfaces, 
but  its  first  partial  derivatives  are  not. 

The  quadratic  fit  for  this  element  requires  six  degrees  of 
freedom . 

(2-12) 

i  - 1 

Where  the  are  basis  functions  and  the  '-fk  are  now  the 
value  of  the  flux  at  corner  nodes  and  the  first  derivative 
with  respect  to  'y  .  Such  that 


Figure  2-2 

Quadratic  Fit  Using  Derivatives  at 
Corner  Nodes 
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node  2  is  at  (0,1,0)  and  node  3  at  (0,0,1)  . 

Over  the  area  of  a  triangle  the  integral  of  natural 
coordinate  powers  is  given  (3:145)  by 
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where  A  is  the  area  of  the  triangie  and  is  equal  to 
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(  2-8) 

The  simplest  interpolant  for  a  triangle  is  the  linear  fit. 
It  express  the  field  variable,  (fi  ,(in  this  case  particle 
angular  flux)  across  the  triangle  as  a  linear  combination  of 
the  flux  at  the  corner  nodes  such  that 


<p  z<f>Uu)  =  <p(*.  4^)  =  £  m 


(2-9) 
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The  partial  derivatives  of  the  flux  are 

r  4  h*  -  k  3.- 

» *■  t  ~ i  it  *  *  -  / 


(2-10) 


t‘a«  ~  "  ~  ~7,  Z’  (2-11) 

where  O’  and  -T*/  are  the  partial  derivatives  of  the  three 

natural  coordinates  with  respect  to  ^  and  M.  respectively. 

Within  an  element  they  are  constant,  but  are  different  for 

each  distinct  triangle  geometry.  Note  that  o<QL  -  —  o 

a  *  a 

It  is  clear  from  this  expression  that  the  equation  is 


satisfied  at  corner  nodes,  that  is  that  3)  =  A.  ,  at  node  1 


element  (3:140).  In  two  dimensions  this  system  is  often 
referred  to  as  area  coordinates,  since  it  can  easily  be  shown 
that  the  natural  coordinates  represent  ratios  of  area. 

Over  a  triangle  one  expresses  the  (  x,u  )  coordinates  in 
terms  of  three  variables  J ^  and  ^  such  that  the  sum  of 

the  natural  coordinates  is  always  one. 


J?,  +  4.  +  =  1 


(2-1) 

The  linear  relation  below  exists  between  the  cartesian  and 
natural  coordinates; 


J,  %  1*  A  V-l  +  A  ^3  “  V 


£,  <u,  +  A  i-  A  ^3  =  M 


Written  in  matrix  form  the  above  relations  become 
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and 


It  is  easily  shown  from  this  expression  that 
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(2-5) 


(2-6) 


2.  of  ^ 

and  that  the  indices  permutate  cyclically.  Note  that  the 

irdinates  of  node  1  in  figure  2-1  are  (  )  -  (1,0,0,) 
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compatibility  requirement  of  elements.  Without  it,  "gaps"  may 
arise  between  elements  from  discontinuous  derivatives,  and 
unsolicited  contributions  to  field  variables  can  arise. 
Completeness  is  the  term  associated  with  the  second 
requirement.  It  insures  that  the  variational  integral  is  well 
defined . 

Standard  finite  element  nomenclature  is  to  express  element 

continuity  as  a  function  of  the  highest  order  derivative  held 

0 

continuous  on  boundaries.  C  continuity  implies  field 

1 

variable  values  are  continuous  on  element  edges,  C 
continuity  has  variable  and  first  derivative  inter-element 
continuity,  an  so  on. 


* 


B.  The  Triangle  and  Two  Dimensional  Interpolating  Functions 
The  two  dimensional  element  chosen  for  this  study  was  the 
triangle.  It  is  a  simple  element  to  refine,  can  be  maneuvered 
easily  to  snugly  fit  irregular  boundaries,  and  can  be 
expediently  described  in  terms  of  its  natural  coordinates 


Figure  2-1 

Cartesian  and  Triangular  Coordinates 
A  natural  coordinate  system  is  a  local  coordinate  system 
that  relies  upon  the  element  geometry  for  its  definition,  and 
whose  coordinates  range  in  value  from  zero  to  one  within  an 
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2 .  Elements  and  Interpo  latlng  Functions 

The  study  of  interpolating  functions,  and  the  elements 
composing  finite  element  meshes  is  an  important  one.  The 
wrong  choices  can  cause  excessive  computations,  or  worse 
prevent  convergence  from  occuring.  The  elements  and 
interpolating  functions  presented  here  are  by  no  means 
inclusive.  Finite  element  literature  on  the  topic  is 
extensive . 

In  this  section  general  requirements  for  monotonic 
convergence  of  the  finite  element  method  will  be  presented 
Then,  the  two  and  three  dimensional  elements  chosen  for  this 
study  are  described.  Finally  the  interpolating  functions  used 
for  each  element  are  given,  and  their  derivations  are 
explained . 


A.  Compatibility  and  Completeness 

Convergence  of  the  finite  element  method  is  guaranteed  if 
the  elements  composing  the  mesh,  and  the  selected 
interpolation  function  satisfy  two  requirements  (3:79). 

Namely,  1)  along  element  boundaries  the  field  variable  and  any 
of  its  partial  derivatives  up  to  one  order  less  than  the 
highest  order  derivative  appearing  in  the  variational 
functional  must  be  continuous  and  2)  the  field  variable,  and 
its  partial  derivatives  up  to  the  highest  order  appearing  in 
the  functional  should  be  represented  in  the  interpolation 
function  as  the  limit  of  element  size  approaches  zero. 

The  first  of  these  requirements  is  the  so  called 
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variable  value  at  that  node  must  be  the  sane  from  each  element 
possesing  that  node. 

5)  Solve  the  System  Equations.  The  resulting  coefficient 
matrix  (referred  to  as  the  global  matrix)  is  symmetric,  banded 
and  positive  definite.  System  of  equations  with  coefficient 
matrices  of  this  type  are  best  solved  by  Cholesky 
decomposition  (7:13)  and  it  is  the  algorithim  used  in  this 
study.  Solution  of  the  system  equations  yields  nodal  values  of 
the  field  variable,  which  together  with  the  interpolating 
function  defines  piecewise  approximations  across  the  domain 
under  scrutiny. 

6)  Hake  additional  Computations  if  desired.  With  respect  to 
the  transport  equation  this  step  is  not  required,  except  in 
the  case  where  penalties  are  desired  for.  automatic  mesh 
refinement . 

0)  Summary 

Operating  on  the  transport  operator  with  the  adjoint 
transport  operator  produces  a  self  adjoint  transport  equation. 
Thi3  equation  can  be  expressed  var ia t Iona  1 ly  as  a  quadratic 
functional,  that  when  minimized  is  equivilent  to  solving  the 
SATE.  The  finite  element  method  is  best  suited  for  solving 
this  type  of  problem.  It  is  a  numerical  technique  that 
approximates  field  variable  values  with  separate  analytical 
expressions,  of  the  same  order,  across  a  mesh  of  small 
interconnected  elements.  The  resulting  set  of  linear  equations 
is  positive  definite,  and  can  be  solved  quickly  by  direct 


means  . 


in  step  3  below.  Additionally,  triangles  are  easily  generated 
and  refined,  and  fit  irregular  boundaries  snugly.  In  general 
one  should  start  with  a  sparse  mesh  composed  of  few  elements, 
this  allows  a  solution  to  be  calculated,  and  regions  of  high 
penalty  identified  for  mesh  refinement.  Constructing  large 
meshes  by  hand  is  a  time  consuming  and  error  prone  process. 
References  for  automatic  mesh  generation  are  listed  in  Heubner 
(3:511) . 

2)  Select  interpolation  functions.  Chapter  two  is  dedicated 
to  this  topic.  Interpolating  nodes  must  be  assigned  to  each 
element,  and  interpolants  chosen.  The  form  of  these  functions 
depend  upon  the  number  of  geometric  nodes  an  element 
possesses,  the  number  of  unknowns  at  each  node,  and  certain 
continuity  requirements  to  be  discussed  in  section  2-a. 

3)  Find  the  element  properties.  The  interpolation  functions 
are  substituted  for  the  field  variable  in  the  functional,  and 
the  integral  is  evaluated.  When  the  resultant  expressions  are 
minimized  with  respect  to  nodal  values,  the  remaining  matrix 
equation  describes  element  properties  in  terms  of  nodal 
variables.  Element  properties  are  expressed  in  terms  of  the 
coefficient  matrices  of  this  equation,  referred  to  in  this 
report  as  the  local  and  nonlocal  matrices. 

4)  Assemble  element  properties  to  obtain  system  equations. 
Uhen  local  and  nonlocal  matrices  are  computed  for  each 
element,  they  are  assembled  globally  to  provide  a  set  of 
simultaneous  linear  equations  with  nodal  values  as  the 
unknowns.  The  foundation  for  this  procedure  lies  on  the  fact 
that  if  a  node  is  shared  by  more  than  one  element  the  field 


than  a  fiaite  difference  rectangular  grid.  Ac  additional 
advantage  of  the  method  is  that  mesh  refinement  can  occur 
easily,  and  that  there  is  a  built-in  indicator  to  dictate 
where  mesh  refinement  should  occur.  Local  mesh  refinement, 
cumbersome  in  a  finite  difference  grid,  consists  only  of 
subdividing  a  finite  element  into  smaller  ones.  Elements  over 
which  this  is  necessary  are  discovered  by  evaluating  the  so- 
called  penalty  function.  Since  the  variational  Integral  is 
minimized,  its  value  over  a  particular  element  is  referred  to 
as  that  element's  penalty.  Elements  where  large  penalties 
occur  are  those  where  the  interpolation  function  fit  is 
poorest,  and  mesh  refinement  is  required.  These  elements  can 
be  subdivided  until  penalties  are  equal  across  the  mesh  and 
further  refinement  fails  to  produce  significant  total  penalty 
reduction. 

Solution  steps 

These  six  steps  are  given  by  Heubner  (3:7)  as  an  orderly 
method  for  obtaining  a  finite  element  solution.  They  are 
described  in  general  terms  below,  and  are  elaborated  upon 
specifically  with  respect  to  the  SATE  in  later  sections. 

1)  Discretize  the  Continuum.  The  first  step  is  to  divide  the 
domain  under  consideration  into  a  set  of  interconnected 
elements.  A  multiformity  of  elements  may  be  used.  The  triangle 
is  very  well  suited  for  two  dimensional  problems,  and  it  is 
the  element  used  in  this  study.  The  ability  to  express 
interpolating  functions  in  terms  of  triangular  natural 
coordinates  simplifies  the  evaluation  of  integrals  appearing 


Across  each  of  these  elements  an  assumed  solution  is 
prescribed,  called  an  Interpolating,  or  approximating 
function.  This  solution  Is  written  as  a  function  of  field 
variable  values,  and  sometimes  the  variable's  spatial 
derivatives,  at  element  nodes.  Solution  requires  choosing 
these  nodal  variables  so  as  to  minimize  the  variational 
statement  of  the  problem  and  satisfy  boundary  conditions.  If 
the  operator  acting  upon  the  field  variable  is  self-adjoint, 
then  equations  describing  the  variable  values  at  element  nodes 
can  be  assembled  globally,  (over  the  entire  material)  and  an 
approximate  solution  to  the  partial  differential  equation  can 
be  calculated  by  solving  the  resulting  set  of  simultaneous 
linear  equations. 

Consider  a  comparison  of  the  finite  element  method  with 
the  well  known  finite  difference  method.  The  finite  difference 
approximation  is  that  a  derivative  can  be  approximated  by  a 
difference  operation  over  a  very  small  interval.  The  resulting 
solution  is  a  set  of  grid  points  at  which  the  field  variable 
values  are  defined.  Finite  elements  assumes  an  analytical 
expression  for  variable  values  over  a  small  element.  This 
approach  results  in  a  piecewise  approximation,  with  field 
variable  values  given  by  separate  analytic  expressions  across 
an  array  of  small,  interconnected  elements,  as  well  as  at 
interpolation  nodes. 

Because  the  finite  element  mesh  is  composed  of  elements, 
they  can  be  put  together  in  a  variety  of  ways.  This  makes  the 
method  well  suited  for  problems  with  complex  geometries.  Small 
elements  can  be  made  to  fit  an  irregular  boundary  much  easier 
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To  be  sure  this  is  correct,  the  expressions  given  by 
equations  (1-9)  and  (1-10)  should  be  satisfied.  Specifically 


-  _1_  1  (*4)X  -  a 

^4  ><4>x 


(1-17) 


should  reduce  to  the  self  adjoint  transport  equation,  and 


^  ‘-Pk 

is  the  condition  required  along  the  boundary.  Straightforward 


(1-18) 


substitution  verifies  that  (1-17)  is  the  SATE,  and  that  the 
natural  boundary  condition  is  the  transport  equation  itself, 
certainly  an  acceptable  requirement. 


Up  to  this  point  the  transport  equation  has  been  recast 
into  a  variational  fora,  and  it  has  been  shown  that  minimizing 
this  functional  is  equivalent  to  solving  the  SATE.  In  the  next 
section  is  presented  background  on  a  numerical  technique  which 
has  achieved  the  most  success  in  solving  this  type  of  problem. 


C)  The  Finite  Element  Method 

The  finite  element  method  is  a  numerical  technique  used  to 
solve  partial  differential  equations.  The  region  under 
scrutiny  is  divided  up  into  a  finite  number  of  elements. 


yields 


(1-11) 


i5^  (V4>  -s)  -sj  1  b 

xr  (^$^4  -JSofo  db  d-12) 

Using  the  definition  of  adjointness  this  functional  can  be 
expressed  as 

X  =  t  £  -  2^3  +*’L)±^  (i-13) 

Setting  the  variation  equal  to  zero,  using  the  definition  of 
adjointness  and  recalling  that ^  of  is  self  adjoint  gives 

«x  =  i  f6  -J.z's)  £*£  £4}  <4t> 

£c-  f  fcp(s^4  -X*s)d0  -o 

=  £4  *x**-*)  ^  =°  (1.14) 

Solution  of  which  is  identical  to  solving  the  SATE. 


The  One  Speed,  One  Dimensional  Functional 

The  one  dimensional,  one  speed,  time  independent  transport 
equation  is  given  (1:76)  by 

/ 

(1-13) 


*„<$  =  €i  fdu' ) 
^  «•  / 


^  ^4 

Tx 

in  the  case  of  isotropic  scatter,  with  no  sources,  where 


jU  ■  cosine(angle  particle  is  traveling) 
p  -  $Cx,  u)  ,  angular  flux 


scattering  cross  section 


In  this  case  ^  -  Z/  \  f  <  -  £5  f  ' 

d-A  ^  T 

This  is  the  form  of  the  transport  equation  chosen  to  test 


the  finite  element  solution  of  the  quadratic  transport 
functional.  The  expression  requiring  minimization  now  becomes 


where  Q 

The  minimum  of  this  functional  is  found  in  an  analogous 
manner  to  finding  the  minimum  of  a  function  in  ordinary 
calculus,  by  setting  the  variation  to  zero. 


.4,  *2. 


ii--f  c  ng £<p  s<£\  ^  j 

Integrating  the  second  term  by  parts 

_fou^v  -  'iruu  —  u; 


o(.  ~  O 


gives 


Cu  =  >  g 

du>  -  ^ d K 
^r<£x 


ns--  S4> 
=S<pK  tL  x 


/^Vy^  .  K'i 

^r=j  /rig  -2.j£l s<pdrUu  + 


(1-7) 


(1-8) 


Since  £"(fi  is  an  arbitrary  admissible  variable  (H:55l)  £"£.  can 

equal  zero  only  if 


*4> 

d-x. 

If. 

|Xl  = 

'*/ 

(1-9) 


JLH  =  0  (1-10) 

The  first  term  above  is  a  simplified  version  of  the  Euler- 
Lagrange  equation  (3:551)  and  is  the  differential  equation 
satisfied  when  is  minimized.  The  second  expression  is 

referred  to  as  a  natural  boundary  condition,  since  it 
specifies  the  solution  form  on  the  boundary,  and  since  the 
functional  can  only  be  minimized  when  it  is  satisfied. 


The  Transport  Functional 

Expansion  of  the  quadratic  functional  (2:15) 
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repeated  here  for  the  purpose  of  document  continuity. 
The  Seif-Adjoint  Transport  Equation  (SATE) 

Adjointness  is  defined  (5:10)  by  the  property 


f  <P^YJr>  =  f  Y  <Lt> 


(  1-4) 


^  T 

where  ^  is  the  adjoint  of  the  operator  .  IfeZ*  *  then^ 

is  said  to  be  self-adjoint. 

Consider  the  operator  /  aC  where  la  the  transport 
operator  and  ^  is  its  adjoint.  Since 

=  =  f Z4>J.b  =f^L<tUC> 

/.  is  self  adjoint.  If  ^  ^  i a  allowed  to  operate  on  the 


transport  equation,  the  resultant  expression 


(1-5) 


is  self  adjoint.  Solutions  of  this  equation  must  satisfy  the 
transport  equation  (2:16). 

/_  (p  —  —  O  is  a  self  adjoint  operator  equation, 

the  solution  of  which  always  satisfies  the  transport  equation. 
This  expression  is  referred  to  as  the  self  adjoint  transport 
equation  ( SATE) . 


B)  Variational  Minimization  of  a  Functional 

The  next  task  described  in  Goff's  thesis  was  to  find  a 
variational  functional,  minimization  of  which  would  be 
equivalent  to  solving  the  SATE.  Before  reproducing  that 
effort,  consider  the  task  of  minimizing  a  functional ,  , 


I-  O,  4>x ;  <L<  d^ 


(1-6) 


M.,  <, 
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1 . Introduction 


A)  The  Transport  Equation 

The  BoLtzman  transport  equation,  written  in  its  general 
in t egr o - d i f f e r en t ia 1  form  is 
v  Yt  - 


S  (isc  <p (*■, 


(1-1) 


where 


%T  is  particle  velocity 
( p  is  particle  angular  flux 

A 

_5"2,  is  particle  direction 

is  the  transport  cross  section  ^  (~r  S7_  tz  *"') 

£  is  particle  sources  S  C  ^  4^  t ) 

<£"  is  particle  energy 
_A/a\ 

2.  is  the  scattering  cross  section  from  energy  E 

__  A  /  A 

to  £  and  angleJl  to  JTi,  .  The  transport  equation  can  be  written 


£ 4>  -*  =  o 

where  the  operator  is  clearly 


(1-2) 


(1-3) 


In  this  formulation,  the  operator  is  non-self-adjoint. 
Finite  element  solutions  of  this  equation  have  been  tried  (1  : 
479)  but  without  a  self-adjoint  operator,  variational  extremum 
principles  do  not  exist,  and  the  finite  element  method's  power 
i3  not  achieved. 

Reformulating  the  transport  operator  into  a  self  adjoint 
form  has  been  accomplished  (2:15)  and  its  derivation  is 


1 


The  basis  functions  are  found  by  requiring  that  at  (1,0,0) 


=*  ,  and  *  .  Four  more  identities  are  found  by  similar 

relations  at  nodes  2  and  3  .  If  the  basis  functions  are 
considered  to  be  the  product 


z'-  ,  _ 

m  6> T 


(2-15) 


where 


/yvi 


A 


(2-16) 


is  a  matrix  of  polynomials,  which  together  represent  a 
complete  quadratic  basis,  then 


(f)  -  ZVTT_  CjT 


(2-17) 


and 


6J  ^ 


(2-18) 


then  since 


2  2  1  ^  5  2  3  1  '  A  3  / 

^  53  ■'“A  5i  | 


the  relation  below  must  hold 
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Now  calculations  are  simplified  since  the  inverse  of  only 
one  3x3  matrix  must  be  found  to  invert  the  6x6  interpolating 
function  matrix.  Flux  at  a  point  in  the  triangle  is  found  by 
evaluating  the  matrix  /W  with  the  natural  coordinates  of  the 
point  in  question,  finding  the  derivatives  of  the  natural 
coordinates  in  the  triangle  with  respect  to  ^  ,  and 
calculating  C£T .  Note  that  Gt  is  a  matrix  of  constants  within 
a  triangle,  but  since  the  depend  on  triangle  geometry, 
will  also  be  different  for  each  distinct  geometry. 

An  unexpected  discovery  prevented  utilization  of  the  above 
quadratic  interpolating  function  in  this  study.  It  is 
described  here  only  because  it  may  be  of  interest  to  other 
researchers,  and  because  it  explains  why  the  more  complicated 
cubic  fit  over  a  triangle  eventually  had  to  be  used. 

For  reasons  to  be  discussed  in  chapter  4,  it  would  be 
extremely  inconvenient  to  construct  a  finite  element  mesh  for 
the  transport  equation  without  the  use.  of  right  triangles.  In 
the  case  of  a  right  triangle  the  derivatives  of  natural 
coordinates  with  respect  to  x  (and  u)  are  constrained  so  that 
two  are  of  equal  magnitude  but  opposite  sign,  and  the  third  is 
identically  zero.  In  this  case  3  of  equation  (2-21)  is 
singular,  and  the  right  triangle  is  in  every  instance 
pathological  for  a  quadratic  interpolant  that  uses  3  fluxes 
and  3  derivatives  as  degrees  of  freedom. 

Used  instead  was  another  quadratic  interpolant,  that  uses 
values  of  flux  at  nodes  and  boundary  centers  to  represent  six 
degrees  of  freedom.  In  this  case 


The  cubic  fit  over  a  triangle  requires  that  10  degrees  of 


freedom  be  specified.  Chosen  were  flux  values,  and  both 
partial  derivatives  at  corner  nodes,  as  well  as  the  triangle 
centroid  field  variable  value.  Assigning  numbers  to  the 
degrees  of  freedom  as  per  below  simplifies  notation. 


A. 

- 


Q\  ^IW  ^2  Li  ^3u  Ou  ^ 


<4  L?z  U>u  'J>r  U),.  <4,  ^  ‘4,0 


(2-28) 


The  basis  functions  are  again  given  by  (2-15)  except  now 


/V* 


J/iS, 

A  A  A~\ 


(2-29) 
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T'T-T'T- r*  *7  ».-t-  r-  i  -  (7  w  p-  »«  ■■ 


/*(V  ,•  I 


1  :vi 


'  %  ^4  *4*!, '■'ft '<*•»> 

~  A  A  a  „  *■  -^j  *  5,  4  ■*-  f  i  ^  ^  $2^'  ^1  ___ 

*Lu  '  j^,V,  -V.4-F,  -^.Cf2  4^-^  H  i-ri'^2. 

4^,  42V,  ^4^,  4  4  4  -4  4  (2-31) 

Evaluating  this  modal  at  each  of  the  10  nodes  yields  the 
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(2-32) 


It  is  not  necessary  to  invert  the  10  x  10  matrix  above  if  it 
is  partitioned  into  9  3x3  matrices  and  several  3x1 
vectors  as  below 
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(2-33) 
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if  1  -[>  ^  JStJ 


and  6J;  is  partitioned  similarly  then 


with 


D  =  A 


Oil 

(z 
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i  Q- 
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A— 
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-y= 

Jr 

A. 

p  -  -*">  i  £ 


(2-34) 


£  -  3 


f;  -  _C 


A-  ^ 

i  --  = 

r  =  -al  2  £ 


(2-35) 


The  chore  of  Inverting  a  10  X  10  matrix  la  simplified.  The 

basis  functions,  and  the  degrees  of  freedom  for  this  cubic  fit 

1 

are  specified.  C  continuity  is  achieved  with  this  fit.  The 
only  drawback  is  that  the  basis  functions  depend  upon  triangle 
geometry,  and  therefor  must  be  recomputed  for  each  unique 
triangle  configuration. 


C.  The  Tetrahedron  and  Three  Dimensional  Interpolation 
Func  tions 

In  three  dimensions  the  simplest  element  is  the  four  node 

tetrahedron.  Four  volume  coordinates,  (  L  ,  L  L  L  )  can  be 

1  2  3  4 

used  to  describe  this  element.  A  straight-forward  linear 
relation  exists  between  x,u,u’  co-ordinates  and  a 
tetrahedron's  natural  coordinates: 


-/ 


and  the  are  the  second  coLumn  of  /?■  .  Similarly  and  ^ 

>C/ 


are  the  third  and  fouth  columns  of  A  respectively. 

Integration  of  tetrahedral  coordinates  over  a  volume  is 
conveniently  given  by  (3s  148) 


Gv ^Lhu:  --  c,v  n  * s/ 

J  +  3)! 


c  Pi-%+r+St3) 

where  V  is  Che  element  volume  given  by 


C 3v  - 


r, 

U,  My 


m!  it-L  ^2 


(2-39) 


(2-40) 


To  prescribe  a  Linear  fit  over  this  element  4  degrees  of 
freedom  must  be  specified.  These  can  be  the  values  of  the 
flux  at  corner  nodes  so  that 


(p  =  ^  U 


(2-41) 


A  cubic  requires  twenty  degrees  of  freedom  in  three 
dimensions.  These  can  be  nodal  values  of  the  flux,  and  the 
three  directional  derivatives  at  each  node  as  well  as  face 


centered  values  of  the  fLux,  as  per  figure  2-6,  with  F  as  the 


field  variable.  Nodes  5,6,7,  and  8  are  face  centered  across 
from  nodes  1,  2,  3  and  4  respectively. 


Figure  2-6 

Numbering  of  20  Degrees  of  Freedom  for  a  Cubic 
Fit  on  a  4  Node  Tetrahedron 


The  basis  functions  for  this  fit  are  again  given  by  (2-15) 

with 
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( 2-42) 


is  now  a  20  X  20  matrix  found  by  inverting 
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where 
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(2-45) 


then 
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M  U 

(2-46) 

/M£  =  ' 

Mr 5  =  Mj  ~  1 

(2-47) 
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(2-48) 
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(2-53) 

and 


The  basis  functions  for  this  fit  are  defined.  They  depend 
upon  element  geometry,  and  require  that  4  separate  4X4 
matrices  be  inverted. 


C.  Summary 

This  study  uses  two  elements  to  construct  finite  element 
meshes  They  are  the  triangle  and  the  four  node  tetrahedon. 
Describing  these  elements  in  terms  of  their  natural 
coordinates  is  s traightf oward ,  and  will  be  seen  to  simplify 
later  calculations. 

Four  interpolating  functions,  one  linear,  two  quadratic 

and  one  cubic  were  investigated  over  a  triangle.  The 

quadratic  that  uses  partial  derivatives  as  degrees  of  freedom 

turns  out  to  be  singular  for  right  triangles,  so  a  quadratic 
0 

C  fit  was  substituted.  Two  fits  were  done  on  the 
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tetrahedron,  a  Linear  and  a  cubic.  Any  fit  that  uses  field 
variable  derivatives  as  interpolants  has  geometry  dependent 
basis  functions.  These  increase  the  number  of  calculations 
required  since  they  must  be  found  for  each  distinct  element 
geometry . 


3  .  The  Case  o  £  no  Scatter 

With  scattering  cross  section  of  zero,  the  expression  to  be 
digitized  (1-16)  becomes 

I-  ~  ^  |  +■  —  S^  1  dxju  I/jIL/P  r®  (3-1) 

'•JrX  i  2  J  =rX' 

Only  particle  streaming  and  absorption  is  occuring.  The  first 

integral  of  (3-1)  is  referred  to  as  the  streaming  term,  since  it 
represents  particle  streaming,  the  second  term  is  called  the 
absorbing  term  for  the  analagous  reason,  and  the  third  term  is 
the  boundary  term.  This  last  integral  results  from  the  cross 
product  of  streaming  and  absorbing  terms  and  is  referred  to  as 
the  boundary  contribution  since  without  it  the  natural  boundary 
conditions  which  arise  from  the  integration  by  parts  in  equation 
(1-12)  are  not  satisfied.  These  three  terms  are  referred  to  as 
local  since  they  fit  field  variable  values  limited  to  the 
triangle  under  scrutiny.  In  thi3  section,  a  description  of 
these  terms'  derivation  and  preparation  for  digitization  will 
occur.  Since  the  quadratic  and  cubic  fits  involve  extremely 
long  derivations,  their  results  only  are  presented  in 
appendices.  Lastly,  a  test  case  to  which  an  analytical  solution 
exists  is  described,  and  numerical  results  of  the  various  fits' 
digitization  are  presented. 

A.  The  Local  Terms 

1.  The  Absorbing  Term.  Recalling  (2-17)  the  interpolating 
function  for  (D 
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Since  the  integration  over  X  involves  3  points,  LI  (  local 
integral  )  (of  dimension  10x1)  and  NLI  (of  dimension  1x10)  must 
be  evaluated  at  x»a,  b  and  c,  then 
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(4-13) 


where  NLM  is  the  (  LQ  X  10  )  (k.  x  l)  non  local  matrix  reflecting 

triangle  i  scattering  to  triangle  j.  The  five  triangle  column 
produces  25  such  non  Local  matrices,  all  of  which  must  be 
assembled  globally,  and  must  be  saved  if  triangle  penalties  are 
desired  as  mesh  refinement  indicators. 


In  appendix  G,  subroutine  LCORD  finds  the  natural  coordinates 
of  the  points  (49  for  Weddle's  n-6 )  needed  on  each  triangle  for 


integration,  and  evaluates  /Yn  and  (needed  to  evaluate  the 


second  scattering  integral)  at  each  of  these  points.  ANING 
perfoms  the  angle  integrals,  finding  matrices  L_I  and  NLI  . 
Finally  SPING  calculates  the  space  integral  across  the  width  of 


the  triangle. 


The  subroutines  are  well  documented,  with  a  list 


Simpson's  method  yieLds 
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where 
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Figure  4-3 

SampLe  5  triangLe  column  with  U  upper,  Lower  and 
center  of  triangLe  3  for  X  *  B 

which,  in  terms  of  the  interpolation  functions  is 

j"/  -*A  5  ^  ^  r ^ 4^3 ^  c/^3 

4?y  if/ 


(4-9) 


=  /SdLTamc 


(4- LO) 


where  )  represents  /n?  evaluated  at  the  natural 

coordinates  of  (  Cy  )  referred  to  as  (  ^  ^  )  ,  and  NLIx  = 

is  the  non  Local  LntegraL  of  triangLe  j  at  x  *  c. 


Similarly,  over  the  local  triangle 


triangles,  was  abandoned. 


C.  Numerical  Evaluation  of  the  Scattering  Integrals 

The  first  attempt  to  evaluate  the  scattering  terms  was  to 
calculate  the  integrals  with  numerical  approximations.  Since 
the  cubic  fit  provides  such  high  accurary,  it  was  felt  that  the 
very  good  streaming  and  absorbing  approximations,  with  less 
accurate  scattering,  would  provide  solutions  properly  reflecting 
the  physics  involved  in  a  problem. 

Simpsons  rule  integration,  Weddle's,  and  Weddle's  rule  for 
n  =  6  were  sequentially  tried.  This  section  will  describe 
Simpson's  integration  for  one  of  the  scattering  terms, since  it 
is  the  simplest  to  write  out.  The  second  integral  and  the  other 
two  techniques  are  straightforward  extensions,  and  appendix  H 
contains  subroutines  used  numerically  to  evaluate  both 
scattering  integrals  with  Weddle's  rule  for  n«6 . 

In  the  case  of  integral  A 


' 

J  du  du '  <X  4>  <fi  ' 


(4-6) 


where  X  ■  (  ~  )  , 

2. 


summation  over  the  5  triangle  column 


of  figure  4-3,  is  represented  as 


=  c<  \  dx 


J 

Maa 


du  ftCx.uJ  ^  |  J,A" 


(4-7) 


whereof  ^1,'  represents  integration  over  u  in  triangle  i 
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problem  to  simplify  evaluation  of  the  integral  in  this  manner. 

If  meshes  are  constrained  columnarly  as  in  figure  4-2,  then 
for  example,  integral  A  is 


o-  _/  _  , 


V.  f 


A 


r 


J 


du  1 


u 


A  A 

'  £  £ 

<■=!  1  -  I 


4:  d' 


(4-5) 


-  /  —  / 


where  n  is  the  number  of  triangles  in  a  column.  Integration  over 
u,  can  proceed  from  the  column's  top  triangle  to  the 


Figure  4-2 

Colummar  Finite  Element  Mesh 

bottom  triangle,  one  element  at  a  time,  halting  at  each  triangle 
to  integrate  over  u'  for  all  elements  in  a  column.  Since  n  is 
the  number  of  triangles  in  a  column,  n*n  angle  integrals  are 
evaluated  per  column.  Each  integral  is  integrated  over  space 
separately  and  results  in  a  non  local  matrix  that  must  be 
assembled  globally.  The  bookkeeping  involved  in  evaluating  the 
two  scattering  integrals  is  simplified.  For  this  reason,  the 
quadratic  fit  using  derivatives  as  finite  element  nodes 
discussed  in  chapter  2,  found  to  be  nonexistent  in  right 
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these  terms  numerically,  and  with  analytical  approximations,  and 
describes  the  results  of  these  efforts.  The  cubic  interpolant  of 
chapter  3  was  used  as  the  finite  element  approximation  for  flux. 


B.  Mesh  Arrangement 

Integration  of  equation  (4-3)  is  over  both  angular  variables 
u  and  u'.  This  proves  to  be  cumbersome.  Consider  integration 
from  x*a  to  x»b  of  Figure  4-1 


Figure  4-1 

An  Unrestricted  Finite  Element  Mesh 
Over  the  Benchmark  Solution  Domain  of 
Figure  3-1 


dx 


(4-4) 


and  d)  are  given  by  differing  cubic  in  t  e  r  po  la  ta  t  ion 
functions  in  each  of  the  7  triangles  composing  this  area. 
Proper  bookkeeping  becomes  a  serious  challenge. 

To  avoid  this,  triangles  are  restricted  to  columns.  While 
this  makes  successive  mesh  refinement  cumbersome,  and  will 
probably  be  awkward  if  the  method  is  ever  extended  to  another 
dimension,  it  is  certainly  appropriate  for  an  early  cut  at  the 


4.  The  Case  of  Isotropic  Scatter 


A)  The  Non  Local  Terms 

When  Scattering  in  allowed  to  occur  (1-16)  must  be  evaluated 
in  its  entirety.  This  equation  can  be  rewritten  as 


tr  \  a  *  c  'A 

Z  J  J 


+  <pj  '+  6  (du'S'ts  (  f.  (4-1) 


2  J 


If  one  neglects  the  three  local  integrals  the  non  local,  or 
scattering  integrals  are  left, 


(4-2) 


V* 


Since  their  is  no  u  dependence  in  the  first  integral,  it  may 
be  integrated  out,  then  a  change  of  variables  from  u"  to  u  may 
occur,  allowing  both  integrals  to  be  combined,  resulting  in 

These  are  the  two  scattering  integrals.  Since  they  Involve 
integration  over  u  and  u',  they  result  in  non-local  terms.  Only 
one  of  them  (  £>  d; 7  )  results  in  a  naturally  symmetric  non-  local 
matrix  after  global  assembelage.  These  integrals  are  referred  to 
as  A  and  B  in  this  report  and  the  code  of  appendix  A.  This 
section  describes  the  preparation  for  digitization  of 


mesh  4  is  worse  than  meshes  2  and  3.  Scrutiny  of  element 


pena 1 ti es 

reveals  that 

thi3  is  due  to 

elements  of  mesh  4  being 

refined  at 

the  area  where  the  largest 

derivative  discontinuity 

occurs  (  x 

*0 ,  elements 

4  and  5  ) .  This 

can  be  considered  as 

Total  Penalty 

Mesh 

linear 

quadratic 

cubic 

1 

.02556 

.00337 

.000398 

2 

.01074 

.000559 

.0000474 

3 

.01433 

.000637 

.0000526 

4  • 

.00716 

.000373 

.0000502 

5 

.00136 

.000067 

.0000158 

6 

.00101 

TabLe  3-3 

Mesh  Penalties  as  Convergence  Indicators 


further  evidence  of  the  cubic  fits'  power,  it  is  flagging  to  the 

0 

programmers  attention  the  nonphysical  boundary  condition;  the  C 
fits  are  not  sophisticated  enough  to  display  the  anomaly. 

E.  Summary 

The  three  terms  which  comprise  the  transport  functional  in 

the  case  of  no  scatter  are  all  local.  Their  derivations  are 

straightforward  and  nearly  trivial  in  the  linear  case.  For 

higher  order  interpolants  the  derivation  is  still  easy  to 

follow,  but  very  long.  When  assembled  globally  these  terms 

represent  the  transport  functional.  Setting  the  variation  of 

this  functional  to  zero  leaves  a  positive  definite  set  of 

simultaneous  linear  equations  that  is  solved  to  find  nodal 

1 

values.  The  cubic  interpolant  (with  C  continuity)  is  more 

0 

powerful,  and  may  be  faster  than  the  C  interpolants,  which 


require  excessive  mesh  refinement  before  converging.  Penalties 
are  powerful  Indicators  of  convergence. 


Two  pieces  of  data  appear  as  anomalies  in  Table  3-2.  The  first 
of  these  is  that  the  cubic  fit  for  mesh  1  appears  better  than 
mesh  2  or  3.  Table  3-3  lists  total  penalties  of  the  meshes,  and 
indicates  that  since  mesh  2  and  mesh  3  have  lower  penalties,  the 
finite  element  fit  is  actually  better  in  these  meshes.  Mesh  1 
was  "lucky"  for  the  cubic  in  that  nodal  values  came  out  so  close 
to  the  analytical. 

For  u>0 

#  of  Triangles  Avg  Perc  Diff  of  Analytic  to  Numeric 


Mesh  /  Mean  Free  Path  Linear  Quadratic  Cubic 


1 

1.33 

46 . 6 

(1.5) 

13.73 

(4.1) 

1.37 

(9.2) 

2 

2.00 

28.4 

(2.2) 

5.55 

(5.5) 

2.83 

(12.1) 

3 

2.00 

29.5 

(2.3) 

5.64 

(6.7) 

2.76 

(12.1) 

4 

3.67 

27.2 

(2.9) 

4.24 

(8.5) 

.84 

(18.8) 

5 

8.50 

8 . 7 

(5.0) 

5.45 

(17.4) 

.55 

(35.2) 

6 

21.33 

3.3 

(14.3) 

e  e  • 

•  •  e 

e  #  e 

e  e  e 

Table  3-2 

Comparison  of  Mesh  Refinement  Required  For  Convergence 
of  Local  Terms.  Average  Percent  Difference  is  From  Analytic 
to  Numerical  Solution  for  u  >  0.  Values  in  Parenthesis  are  CPU 
Seconds  of  Runtime  on  a  Vax  11-780,  Unix  Berkely  1.2,  During 
Periods  of  Moderate  to  Almost  Heavy  Use. 


Likewise  mesh  5  for  the  quadratic  fit  appears  worse  than  less 
refined  meshes.  The  total  penalty  bears  out  that  mesh  5  is  a 
better  fit.  These  and  other  similar  experiences  emphasis  two 
points  that  should  not  be  neglected  with  the  finite  element 
method.  The  first  of  these  is  that  element  penalties  can  be  as 
good  a  measure  of  convergence,  or  better,  than  comparing  nodal 
values  to  some  "exact"  solution.  Secondly,  the  meshes  used  in 
this  3tudy  are  not  necessarily  successively  refined.  Without 
this  type  of  refinement,  steady  convergence  of  values  to 

the  exact  solution  may  not  be  observed  (3:79). 

One  anomaly  appears  in  table  3-3.  That  is  the  cubic  fit  for 


refinement.  The  difference  is  clearly  caused  because  flux 


U»1  X 

Analytic 

Linear 

Quadratic 

Cubic 

.375 

.  6873 

.  6360 

.  6904 

.  6865 

1 . 000 

.4724 

• 

.  3834 

.4282 

.4707 

1 . 500 

.  2231 

.1405 

.1707 

.  2230 

3.000 

.  0498 

.  0270 

.0302 

.0512 

CPU  units 

•  •  • 

2.9 

8.5 

18.8 

sec's) 

Table  3-1 

Comparison 

for  Mesh  4  of  Accuracy 

and  Run 

Times 

for  3 

Interpolation 

Fits  on  a 

Vax  11-780 

( unix , 

Berkely  1-2) 

for  u*l 

1  0 

spatial  derivatives  are  held  continuous  in  the  C  fit.  The  C 
quadratic  fit  is  only  slightly  better  than  the  linear.  Consider 
the  expansion  of  1-17  and  1-18  in  the  no  scattering  case. 


-2^  ^4> 


-  o 


( l-17a) 


ft* 


-  O 


( 1- 18a) 


These  are  the  differential  equations  being  satisfied  when 
the  variational  functional  is  minimized,  and  the  natural 


boundary  condition.  Three  quantities  in  these  equations  must  be 

2  0 

approximated  by  finite  elements,  §  ,  \ £  and  ^  41  .  C 

b  *  t*2- 


continuity  holds  only  one  of  these  continuous  across  elements 

1 

boundaries.  The  C  fit  holds  two  of  the  three  continuous  and  as 


a  result  converges  faster.  This  analysis  further  indicates  that 
2  3 

a  C  fit  would  converge  even  faster,  and  a  C  fit  would  be  no 
2 

better  than  a  C  ,  since  no  higher  order  terms  are  needed  to 


satisfy  these  equations. 


Q  -  AO  -2  ~sm"  yU  >  O 


and 


0  =  0 


^  o 


(3-19) 

(3-20) 


Since  most  streaming  is  occuring  near  u  *  1,  meshes  were 
refined  more  in  that  area.  Also,  notice  that  the  boundary- 
conditions  are  not  physical,  there  is  a  discontinuity  of 
derivatives  along  the  u*0  line.  For  now  one  must  realize  that 
this  discontinuity  is  a  source  of  error  that  becomes  apparent 
for  u  near  zero. 


C.  Results 

The  area  of  the  Benchmark  problem  was  discretized  with  6 
separate  meshes  for  the  computer  runs.  All  are  drawn  and  listed 
in  appendix  F.  Hashes  1  through  5  are  those  used  by  Goff 
(2:114)  while  mesh  6  is  a  very  well  refined  mesh  of  80  triangles 
in  3  mean  free  paths.  Mesh  2  and  3  consist  of  the  same  number  of 
triangles,  but  with  a  different  pattern,  to  test  sensitivity  to 
element  orientation.  Table  3-1  shows  a  comparision  of  interpola¬ 
tion  fit  accuracy  with  the  analytical  solution  for  mesh  4. 

Table  3-2  Lists  the  average  nodal  percent  difference  from  the 
analytical  solution  to  allow  a  comparison  of  the  degree  of  mesh 
refinement  required  for  convergence.  It  is  seen  from  these 
results  that  the  cubic  fit  is  extremely  powerful.  Run 
times  should  not  be  taken  as  absolute,  but  it  appears  that  the 
price  paid  in  terms  of  extra  calculations  for  the  more  accurate 
fit  is  not  extreme.  The  linear  fit  converges  as  finite  element 


theory  says  it  will,  but  only  with  an  excessive  amount  of  mesh 


performed  for  a  triangle,  the  resultant  quadratic  form  can  be 


written  as 


-  — ] 

M  fir  -f.  /U  0  -I- /US  (JT  ~  jf. 


bT  ML  6T  ^ 


(3-17) 


where  ML  is  the  local  matrix.  The  global  assembelage  of  these 
local  matrices  results  in  a  quadratic  form  for  the  variational 
integral  over  the  entire  area  under  scrutiny.  When  minimized, 
the  remaining  set  of  simultaneous  linear  equations  has  a 
coefficient  matrix  that  is  always  positive  definite,  and  is 
solved  by  cholesky  decomposition  in  this  study. 


B.  The  Test  Case 

The  problem  chosen  to  test  the  digitization  of  the  local 
terms  was  a  monoenerge tic  lambertian  (flux  ■  cosine  of  the  angle 
it  is  traveling)  source  of  particles  incident  upon  an  absorbing 
only  slab  3  mean  free  paths  thick.  The  slab  is  surrounded  on 
both  sides  by  a  vacuum  so  there  is  no  returning  flux  from  the 
right  boundary. 


M 


Figure  3-1 

Benchmark  Problem  Description 

In  this  case  the  transport  equation  is 

M  IS  ^  ~  O  (3-18) 

is  < 


The  solution  is 


interpolant  evaluation  are  in  appendix  C  . 


3.  The  Streaming  term. 

Expansion  of  M.  yields  six  integrals  which  must  be  evaluated, 
JJ.  ~  At,  4.  +  +~  ^3 

■h  4^3 


-i-  p  A  j^y’  '  5  4T  +pA  4*4' 

^ jd 4  AjV/a*^  2^  +  pA  2A4,Mr^/7_  M\y  ^ 

+p*  2*a,uz  tjz  *2?^  +  P4  V3 


(3-11) 


— >y  -/y 


(3-12) 


=  4^4:  £*4*  +/^12  4-/^3  6r  <£ 


•=  a  ^  t=J  ^ 


(3-13) 


(3-14) 


When  linear  interpolants  are  substituted  the  streaming  term 


l 

TJ.1  5,^  9i 


£?  =  — *  e°i  3i  91,3 


).9j  33i 


f  z.  /J,(uxrM -u)  +-  +  MiC^Z  tr/U>)  J 

Derivation  of  this  term  for  the  quadratic  and  cubic  fits  are 
the  most  lengthy  of  all;  results  of  this  effort  are  in  appendix 


(3-15) 


(3-16) 


4.  The  Local  Matrix.  The  sum  of  these  three  terms,  evaluated 
over  a  triangle  in  the  mesh,  is  the  value  r c  the  variational 
integral  over  that  element.  After  these  calculations  are 


3 


N  H 


The  Boundary  term  can  be  written  as 


—  -  £  2yU  °p-  is  =  c  /?  -<w  O'  6  /vh  ^  ^  A,  7“ 


(3-6) 


using  (2-17)  and  since 


-  ‘±  kS  /w  ^ 


(3-7) 


the  Boundary  term  can  be  expressed  as  the  sum  of  three  integrals 


_  S\~- 

—  2_  —  o'  A  //.'  2  '  V  T 


/  6  r  ^  ^  (gT* 


(3-8) 


in  the  linear  case,  evaluation  of  the  integral  yields 
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(3-9) 


The  resultant  matrix  will  not  be  symmetric,  but  since  for 
each  quadratic  form,  only  one  symmetric  matrix  exists  (4:342), 
the  boundary  term  can  be  symmetrized.  That  is 


This  term  is  much  more  complicated  to  encode  than  the 
absorbing  term,  since  it  involves  derivatives  of  natural 
coordinates  with  respect  to  x,  which  are  different  for  each 
separate  geometry.  The  results  of  quadratic  and  cubic 
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(3-2) 


The  absorbing  term  can  be  written  as 


r  r 

-  1  /  C  '-A 

"2.  . 


(p  r  id/3  (jJ  /W  a£  d?T~ 


r, 


^•fr  6_J"  p./f  W/w  !6i  if 

-2.  "■  ^  "  — 
L_  _ 


Since  (j?  and  6J"  are  constant  within  an  element,  they  can  be 
removed  from  the  integraL.  Evaluation  of  the  term  then  reduces 
to  taking  the  product  of£*/w  and  analytically  performing  the 
integral . 


In  the 

linear  case  6 T  ■ 

1  X  and 

‘i,  " 

“jX 
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which , 

using  (2-7)  is 
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where  A 

is 

the  area  of 

the 

triangle 

i . 

The  absorbing  term  can 

now 

be 

written  as 
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(3-3) 


(3-4) 


(3-5) 


It  is  precisely  the  above  expression  in  brackets  which  is 
digitized  and  evaluated  for  each  element. 


0  1 

With  the  C  quadratic  fit,  and  C  cubic  previously 
described  in  chapter  2,  the  term  derivation  is  analagous,  expect 
tha t  A-  in  these  cases  is  of  order  6  and  10  respectively. 
Appendix  B  has  results  of  these  computations. 

2.  The  Boundary  term 


of  variables  included  in  the  appendix. 


Note  that  LI  and  NLI  are  actually  the  same  integral,  just 
over  different  triangles. 

J'.jw.  4  =  ? 

J1”"  u 1  dp  —  M  4J*  ^  ^ 


Subroutine  ANING  takes  advantage  of  the  fact  that  L_I  *  NLI 
and  calculates  over  every  triangle  in  the  mesh, 

storing  it  in  memory  to  be  recalled  when  needed.  Simultaneously 


it  calculates 


u.  u. 


,  the  only  other  angle  integral 


needed,  storing  these  in  ILFD , 


c* 


D.  Cubic  Analytical  Approximation 

Since  (p  and  cj)^  are  approximated  with  cubic  functions,  the 
scattering  integrals  of  eqn  (4-3)  ,  which  integrate  the  products 
of  (jp  and  q£>  ^  are  integrating  a  hexadic.  Explained  in  this 
section  is  an  attempt  to  substitute  another  cubic  for  the 
"exact”  sixth  order  fit  required  by  4? '  and  4^  $ 

In  three  dimensions,  the  local  and  non  local  triangles  map 
out  tetrahedrons.  Four  possible  cases  can  occur,  depending  upon 
the  orientation  of  the  local  and  non  local  triangles  as  depicted 
in  Figure  4-4.  Case  2  and  case  4  result  in  pyramids,  which  can 
split  along  their  center  into  two  four  node  tetrahedra  each,  and 
integration  can  then  be  performed  separately  over  each 
tetrahedron,  and  summed. 


Consider  the  first  scattering  integral 

J~  A  K  <kM  <Lm'  OC  (p  <p/ 


(4-14) 

(4-15) 
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F,  the  product 


can  be  approximated  in  three  dimensions  with 


a  cubic  polynomical  of  complete  basis  as 


j—  —  /y/\  T"  -f~ 


(4-16) 


where/y/i  and  GT  are  as  given  in  (2-41)  and  (2-44)  respectively. 

F  can  be  considered  to  be  a  column  matrix  of  twenty  10  X  10 
matrices,  one  representing  F  at  each  of  the  twenty  nodes  of 
figure  2-6.  For  instance,  case  1,  node  1  is  the  intersection  of 
nodes  2  (local)  and  node  1  (primed).  f  is  then 


/  .  / 


-fj r  4)  '  =  C£  /w  fa'  (yr  ld 


(4-17) 


(4-18) 


evaluation  of^  at  (0,1,0)  and  /vh  at  (1,0,0),  and  carrying  out 
the  multiplication  is  guaranteed  to  yield 

~o  cd  o  o  o  o  a  o  o  o 

I  Q  o  CD  o  O  O  o  O  O 

o  o  o  O  O  CD  o  O  O  o 

0  o  a  o  O  o  o  a  o  a 

0  a  CD  0  C  CD  o  O  o  cd 

^<-1  o  o  o  CD  O  O  O  O  CD  CD  (4-18) 

O  O  CD  CD  CD  CCD  CD  3  CD  O 

o  O  CD  o  o  CD  OO  CD  CD 

C  O  CD  O  CD  CD  CD  CD  O  O 

OO  CD  O  CD  CD  CD  D  O  CD ^ 

since  both  and  <^>  /  are  finite  element  interpolation  nodes. 

With  this  cubic  approximation,  the  second  scattering  integral 
is  only  slightly  more  complicated. 


-2,CW« 

if  the  same  cubic  approximation  is  made  for  G»  (p*  {j) 


(4-19) 
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:ri. 


then  four  integrals,  caused  by  the  expansion  of  u,  must  be 
eva 1 ua  ted  . 

The  5  and  ^  *5  are  in  most  instances  trivial,  and  can  be 
written  down  by  inspection  for  each  of  the  four  cases.  This  is 
done  in  appendix  E. 

The  integrations  are  simple  compared  to  those  done  for  the 
streaming  case  since  there  are  no  cross  products  of  and  ^*1^. 

E .  The  Test  Case 

Chosen  to  test  the  numerical  and  analytical  evaluations  of 
the  integrals  was  the  same  domain  as  in  the  no  scattering  case, 
with  the  region  depth  under  scrutiny  varying  from  one  to  five 
mean  free  paths.  Graciously  provided  by  Dr.  Shanlcland  was  a 
spherical  harmonics  solution  of  the  problem  using  up  to  46 
legendre  polynomials.  Results  of  these  calculations  are  listed 
in  appendix  G  for  scattering  cross  sections  corresponding  to  c 
of  .5  and  .9  where  D*.  Shankland  used  as  a  lower  right 

boundary  condition  no  return  flux  at  infinity.  Therefore,  the 
lower  right  boundary  used  in  the  finite  element  solution  is  the 
Pn  angular  flux  for  u  <  0  at  1,2, 3, 4  or  5  mean  free  paths, 
depending  upon  the  depth  of  investigation  desired.  In  this  case, 
there  are  no  sources  in  the  region  under  scrutiny,  and  the 
coupled  Pn  equations  are  solved  with  a  Green's  function. 

The  Lambertian  source,  depicted  in  figure  3-1  is  non 
physical.  The  derivative  discontinuity  at  x*0  is  very  difficult 
to  approximate  with  a  finite  polynomial  series.  The  expected 
solution  for  the  lambertian  at  this  spatial  point  for  the  cases 
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of  c».5  and  c=.9  would  be  similar  to  figure  4-5,  with  more 


backscatter  in  the  c».9  case  than  the  c».5  .  As  a  result,  the 

approximating  function  generated  by  the  legendre  polynomials 
changes  less  rapidly  about  u*0  for  c*.9  than  for  c*.5,  and  can 
be  constructed  with  less  polynomials. 


Figure  4-5 

Expected  Solution  at  x"0  for  the  Lambertian  Source 
and  Expected  Legendre  Polynomial  Approximations 

Used  as  finite  element  left  hand  boundary  conditions  are  the 
legendre  polynomial  values  of  angular  flux  in  appendix  H  for  X*0 
and  u>0  .  With  these  boundary  conditions  the  finite  element 
solution  was  tested,  and  its  results  compared  to  the  spherical 
harmonics  solution  throughout  the  regions,  of  varying  depth,  and 
with  varying  cross  sections,  under  scrutiny. 

D.  Results 

Penalty  functions  are  guaranteed  to  be  positive  with  this 
method  since  the  value  of  the  functional  is  the  integral  of  a 


quantity  squared. 

I =  i  f  ^  (**4  - 0  C^4  ■ s) 


(1-11) 


If  a  negative  penalty  occurs,  it  is  an  indicator  of  error.  In 
the  case  of  no  scatter,  every  penalty  was  greater  than  zero 
because  flux  was  approximated  as  a  cubic,  and  the  integration 
was  exact.  In  the  scattering  case,  numerical  evaluation  of  the 
scattering  integrals  is  hardly  exact,  nor  is  approximating  a 
hexadic  with  a  cubic,  and  then  analytically  performing  the 
integral.  Negative  penalties  could  occur,  and  they  would 
indicate  error  in  evaluating  the  scattering  integral. 

The  global  matrix  is  guaranted  to  be  positive  definite.  In 
the  streaming  case  it  always  was  positive  definite,  again  due  to 
the  exactness  of  the  integration.  In  the  scattering  case,  a  non 
positive  definite  global  matrix  is  another  indicator  of  error  in 
scattering  integral  evaluation.  This  type  of  matrix  could  be 
solved,  and  the  solution  might  be  fairly  accurate,  but  a 
desirable  chara t er i s tic  of  the  finite  element  method  is  that  the 
resulting  set  of  linear  equations  has  a  positive  definite 
coefficient  matrix,  since  it  can  be  solved  quicky  and  accurately 
by  direct  means.  Other  than  positive  definite  matrices  solution 
techniques  must  rely  upon  iterative  solution  methods,  or  very 
long  direct  schemes,  therefore  this  study  insists  that  the  means 
used  to  evaluate  scattering  results  in  positive  definite  global 
coefficient  matrices. 

Numerical  Results 

Both  negative  penalties  and  non  positive  definite  matrices 
were  common  with  the  three  numerical  techniques  used.  Simpsons 
rule  was  used  first.  For  simple  meshes  (meshes  1-4)  positive 
definite  global  matrices  occured.  For  any  further  refinement, 


error  in  the  integration  grew,  and  the  matrix  became  non 
positive  definite.  Since  Weddle's  rule  fits  a  cubic,  it  was 
tried  next,  and  more  refinement  could  occur  until  the  same 
effect  happened.  Weddle's  rule  for  n*6  was  the  last  numerical 
technique  tried,  and  its  error  causes  non  positive  definite 
matrices  with  around  15  triangles  per  mean  free  path  at  c«.5  . 

In  table  4-1  are  results  for  Weddle's  N*6  rule  in  mesh  D  with  a 
depth  of  3  mean  free  paths  and  c«.5.  It  shows  that  the  method 
does  not  provide  acceptable  accuracy.  It  seems  odd  that  mesh 
refinement  would  increase  error.  What  is  occuring  is  that  as  the 
number  of  triangles  is  increased,  the  numerical  integration  must 
be  performed  more  often.  The  error  accumulates  until  it  destroys 
the  global  matrices  positive  definiteness.  This  is  even  more 
apparent  if  c*.9  is  used;  the  scattering  integral's 
contributions  are  greater,  and  even  less  refined  meshes  produce 
negative  definite  matrices.  Numerical  integration  of  the 
scattering  integrals  holds  no  potential.  Orthogonal  relations 
could  be  tried,  they  are  used  successfully  in  the  Sn  method,  but 
they  would  probably  not  be  the  solution.  In  the  Sn  method, 
iteration  throughout  the  mesh  must  occur  to  reach  the  proper 
solution.  The  finite  element  technique,  as  formulated  in  this 
study,  is  not  adaptive  to  iterative,  or  "marching”  methods. 

Cubic  Approximation 

Comparison  of  finite  element  and  PN  solution  for  c*.5  and 
c*.9  were  conducted.  Two  types  of  boundary  conditions  were 
tried.  First,  only  fluxes  were  specified  and  secondly  fluxes 
and  its  derivatives  with  respect  to  u  were  specified.  Derivatives 
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Table  4-1 

Weddle's  Rule  n*6  Results  for  Mesh  D,  with  Range-3.0, 
u  Derivatives  and  Fluxes  Specified  on  the  Boundary,  c».5 
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were  found  by  the  use  of  difference  equations  on  the  Pn  data.  A 
total  of  four  meshes,  A  through  D  were  used  (appendix  F).  A,B, 
and  C  meshes  all  have  a  depth  of  3  mean  free  paths.  Mesh  D  depth 
was  varied  from  1  to  5  mean  free  paths.  Since  data  from  these 
meshes  is  extensive,  selected  output  is  displayed  in  appendix  A, 
and  results  are  summarized  in  table  4-2. 

The  results  of  table  4-2  indicate  that  convergence  is 
occuring  for  c*.9  data  only  after  excessive  mesh  refinement. 

The  data  for  c*.5  indicates  convergence  of  the  finite  element 
code  is  occuring,  but  not  to  the  Pn  solution.  Specifying  u 
derivatives  speeds  up  convergence.  Penalties,  expounded  as 
being  so  important  in  chapter  3,  appear  to  carry  no  useful 
information. 

Both  boundary  conditions  are  appropriately  specified.  It  is 
customary  in  the  widely  used  Pn  and  Sn  transport  codes,  to 
specify  only  fluxes.  In  general  these  codes  do  not  directly  use 
the  u  derivatives  on  boundary.  However  if  the  flux  is  known  as 
a  function  of  u  at  a  specific  spatial  location,  then  certainly 
the  flux  derivative  with  respect  to  u  is  known  at  that  point. 
Since  the  finite  element  code  uses  C as  an  interpolation  node, 
specifying  its  value  on  the  boundary  is  appropriate,  and  can 
only  speed  convergence  to  the  same  solution. 

The  lack  of  exactness  in  scattering  integral  evaluation  has 


destroyed  element,  and  total  penalty  usefulness.  Consider  the 

origin  of  a  particular  element's  penalty,  pen(i) 

b et+gy* 

=  ^  A JLM(  c  j)  '±i- 

1  1  J~vf 

J  aa/i/*r»w  a.  (4-21) 
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Table  4-2 

Cubic  Approximation  of  Scattering  Integral  Results 
Compared  with  Legendre  Polynomial  Solution  for  c».5  and 
c».9.  Average  Percent  Difference  of  Nodal  Values  Other  than 
Those  Specified  as  Boundary  Conditions,  with  flux  only  as  a 
Boundary  Condition.  Values  in  Parenthesis  are  Same  Meshes 
with  Flux  and  u  Derivative  Specified.  *  is  Non  Positive 

Definite  Global  Matrix. 

exactly  this  that  constitutes  penalties  in  the  case  of  no 
scatter.  The  scattering  component  is  expected  to  be  negative, 
as  can  be  seen  by  evaluating  the  signs  of  the  integral 
coefficient  constants  of  equation  (4-3),  ~ 

are  both  Less  than  zero.  The  sum  of  the  streaming  and 
scattering  penalty  contributions  should  remain  positive  however, 


52 


since  they  represent  the  square  of  a  quantity.  When  scattering 
evaluation  is  with  error,  its  penalty  contribution  can  grow  too 
large,  and  an  element's  penalty  drops  below  zero.  If  this 
happens,  summing  of  element  penalties  for  a  total  mesh  penalty 
leads  to  misleading  information.  It  was  thought  that  the 
magnitude  of  a  penalty  might  carry  the  desired  information  on  a 
fit's  correctness,  so  the  sum  of  penalty  absolute  values  was 
computed.  Comparison  of  this  data,  displayed  in  table  4-2  also 
lacks  the  desired  information.  Negative  element  penalties  are  in 
every  instance  associated  with  triangles  where  a  large  amount  of 
scattering,  and  a  small  amount  of  streaming  is  occuring.  Because 
the  scattering  integral  evaluation  is  inexact,  the  penalty 
function  has  lost  its  value. 

Graphs  of  figures  4-6  and  4-7  compare  finite  element 
solutions  with  the  legendre  polynomial  solution,  for  various 
triangle  densities  and  numbers  of  legendre  polynomials  being 
used.  All  finite  element  computations  on  the  graphs  were  done 
with  mesh  D,  varying  the  depth  to  change  triangle  densities.  It 
appears  from  this  data  that  with  more  legendre  polynomials  the 
finite  element  and  Pn  solutions  would  be  exact.  Convergence  is 
faster  in  the  c*.9  case  because  the  larger  backscatter  creates  a 
smoother  flowing  function,  able  to  be  approximated  with  fewer 
legendre  polynomials  than  the  more  rapidly  changing  c=.5 
solution.  Boundary  conditions,  used  in  the  finite  element  code 
as  specified  by  the  spherical  harmonic  solution,  have  not 
settled  down  yet  either,  as  shown  in  table  4-3. 

Based  upon  this  information  it  appears  that  the  finite 
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Table  4-3 

Pn 

Predicted 

Flux  Rate  of 

Change  of  Selected 

Boundary 

Points  Over  Polynomials  30-46 

element  calculations  are  succesfully  predicting  angular  flux. 
Under  these  circumstances  it  is  difficult  to  determine  the 
amount  of  refinement  required  for  convergence,,  but  it  appears 
that  if  fluxes  only  are  specified  on  boundaries  the  method 
converges  with  15  to  20  triangles  per  mean  free  path.  If  fluxes 
k  and  derivatives  are  specified,  convergence  occurs  with  10  -  15 

elements  per  mean  free  path.  Less  angle  refinement  is  also 
required  if  derivatives  are  specified  on  boundaries.  This  is 
approximately  the  same  degree  of  refinement  as  an  S4  calculation 
with  two  spatial  nodes  per  mean  free  path.  Positive  definite 
matrices  are  not  guaranteed  (  as  in  the  case  of  mesh  B.). 

The  only  difference  between  mesh  C  and  D  is  refinement  over 
angle  in  the  first  and  last  columns.  Close  analysis  of  table  4-2 
data  shows  that  this  angle  refinement  is  more  important  than 
spatial  refinement  when  fluxes  only  are  used  as  boundary 
conditions.  This  is  further  indication  of  scattering  term 
inexactness.  The  scattering  calculations  error  can  be  estimated 
from  c».9  data  of  table  4-2.  With  46  triangles  per  mean  free 
path,  the  average  nodal  percent  difference  of  1.117.  can  be 
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Pen  -  element  penalty 

Gl,  G2,  G3  -  derivatives  of  triangular  coordinates  w.r.t.  space 
FI,  F2,  F3  -  derivatives  of  triangular  coordinates  w.r.t.  angle 
X1,X2,X3,U1 ,  U  2  ,  U  3  -  specific  coordinates  of  the  triangle  under 
scrutinies  geometric  nodes 

V  -  Array  storing  the  integral  of  the  twenty  tetrahedral 
coordinate  combinations  which  together  form  a  complete 
basis  for  a  cubic  in  three  dimensions  (2-41).  Row  two  has 
the  integral  of  times  (2-41),  row  three  times  (2- 

41)  ,  rows  4  and  5  contain  and  times  (2-41) 

integrated  over  tetrahedral  volume  respectively. 

SGM  -  M  5 ,  M  6 ,  M  7 ,  and  M8  of  (2-42)  and  Ml  8  of  (2-43) 

E,  F,  G,  -  arrays  of  dimension.  4  storing  the  derivatives  of  the 
four  tetrahedral  coordinates  w.r.t.  space,  incident  angle  and 
scattered  angle  respectively 

H  -  the  basis  functions  for  each  of  the  5  scattering  integrals 
(expansion  of  u  requires  that  the  second  integral  be  done 
four  separate  times) 

SA  -  the  first  scattering  integral  matrix 
SB  -  the  second  scattering  integral  matrix 

Integers  Passed  as  Arguments 
N  -  number  of  nodes 
TRI  -  Local  triangle 
TRIP  -  non  Local  triangle 
STRIA  -  number  of  triangles 
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Appendix  A  -  Program  Listing 

Glossary  of  Variables 

Variables  Passed  as  Common 
MG  -  Global  matrix 
ML  -  local  matrix 
NLM  -  non  local  matrix 

GT  -  matrix  of  interpolating  function  constants 

Variables  Passed  as  Double  Precision  Arguments 
Cordnd  -  cartesian  (x,u)  coordinates  of  finite  element  nodes 
Phi  -  angular  flux 
Areas  -  triangle  areas 
MA  -  absorbing  matrix 

SCI,  SC2  ,  ...  SC6  -  coefficients  of  streaming  matrices  per 

appendix  D 

SRI,  SR2 ,  ...  SR6  -  per  appendix  D,  row  matrices  to  augment  SC 

matrices 

BC1 ,  BC2 ,  BC3 ,  BR1 ,  BR2  ,  BR3  -  coefficients  of  boundary 
matrices  per  appendix  C 
MB  -  boundary  matrix 
MS  -  streaming  matrix 

DRVS  -  matrix  of  derivatives,  overlayed  on  boundary  term 
coeficients  per  appendix  C 

Range  -  depth,  in  mean  free  paths,  of  region  under  scrutiny 
SIGMAT  - 
SIGMAS  -  2.  s 
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Any  interpolant  that  uses  field  variable  derivatives  as 
finite  element  interpolation  nodes  has  basis  functions  that  are 
geometry  dependent,  and  increases  calculations  significantly. 
This  type  of  interpolant  was  accepted  for  the  local  terms 
because  the  increase  in  accuracy  made  up  for  the  extra 
calculations.  It  is  not  necessary  to  use  a  geometry  dependent 
interpolating  function  for  the  nonlocal  terms,  and  this  type  of 
approximation  holds  no  accuracy  benefits.  Exact  scattering  term 
evaluation  is  possible  with  geometry  independent  interpo lan ts , 
and  is  recommended  as  any  subsequent  study's  first  effort. 


interpo lants  . 


The  streaming  results  clearly  show  the  penalty  function's 
usefulness.  Not  only  is  the  global  penalty  a  faultless 
indicator  of  accuracy,  element  penalties  may  be  used  to  dictate 
where  local  refinement  should  occur. 

Two  methods  were  used  to  evaluate  the  scattering  terms,  and 
analysis  of  their  results  leads  to  a  proprosal  for  a  third 
integration  technique,  which  should  be  both  more  efficient  and 
accurate.  Numerical  techniques,  of  accuracy  up  to  Weddle's  for 
n*6  were  unsuccessful.  Cubic  approximation  of  the  hexadic 
scattering  integral  was  accurate,  required  slightly  more 
refinement  than  expected,  and  appears  to  be  computationally 
excessive.  Worst  of  all,  the  inexactness  of  the  scattering 
integral  evaluation  destroys  penalty  value,  and  does  not 
guarantee  positive  definite  global  matrices.  Chapter  5 
describes  proposed  exact  hexadic  integration,  with  geometry 
independent  basis  functions  that  should  significantly  reduce 
computations,  return  penalty  usefulness,  and  insure  positive 
definiteness.  With  exact  scattering  integral  evaluation, 
accuracy  equal  to  the  streaming  case  should  be  achieved  with 
comparable  mesh  refinement. 

Extensions  to  other  than  isotropic  scatter  will  be 
straightforward.  If  the  scattering  kernel  is  expanded  in  terms 
of  a  legendre  polynomial  series  as  it  ordinarily  is,  the 
scattering  integrals  would  be  slightly  more  complicated,  but 
achievable.  Integration  with  dx,  du  and  d.u '  over  a  four  node 
tetrahedron  wouid  still  result,  only  the  form  of  the  function 
being  integrated  would  be  changed. 
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6.  Conclusion 


The  finite  element  method  has  been  very  successful  in  a 

variety  of  fields.  It  was  felt  that  since  the  self  adjoint 

reformulation  of  the  transport  operator  could  be  expressed  as  a 

quadratic  functional,  finite  elements  could  be  applied 

succesfully  to  transport  problems.  Concisely  stated,  the  result 

of  this  study  is  that  the  method  works,  and  that  it  appears  to 

hold  potential  for  very  accurate  solutions  with  moderately 

refined  meshes.  The  present  digitization  of  the  method, 

described  in  this  document,  and  written  in  appendix  A,  bears 

improvement,  both  in  accuracy  and  computational  efficiency. 

It  was  found  that  with  linear  interpolants  the  method 

converged  in  the  case  of  no  scatter,  with  around  25  triangles 

per  mean  free  path  for  u>0.  Linear  interpolants  were  not  tested 

in  the  scattering  case,  but  straightforward  extension  of  the 

streaming  results  suggests  that  at  least  50  triangles  per  mean 

free  path  will  be  required  to  reach  an  accurate  solution.  This 

0 

is  an  enormous  amount  of  refinement.  The  C  quadratic  fit  was 

only  slightly  better.  Unfortunately  columnar  mesh  restriction 

1 

destroyed  a  semi  C  fit,  and  cubic  interpolants  were  used. 

These  are  very  powerful  in  the  streaming  case,  achieving 

accuracy  of  greater  than  9  97.  with  around  4  triangles  per  mean 

free  path.  Codes  used  in  this  study  where  not  written  with  the 

intention  of  comparing  speeds.  Run  time  comparisons  are 

therefore  not  absolute,  but  they  do  indicate  that  the  more 

accurate  fit  is  not  computationally  excessive,  and  may  even 

0 

require  less  cpu  time  to  converge  than  either  of  the  C 
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where  k  and  Ay  in  this  instance  represent  the  dimension  (10) 


distinct  basis 

functions  and 

their  spatial  derivatives  ■ 

cubic  fits  over 
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Table  5-1 

84  Polynomials  for  Three  Dimensional  Hexadic 


C.  Summary 

Results  of  chapter  4  dictate  the  need  for  exact  scattering 
integral  evaluation.  The  hexadic  using  flux  only  as  degrees  of 
freedom  will  integrate  exactly,  and  probably  reduce 
computations.  Time  precluded  digitization  of  this  fit,  and  it 
is  recommended  as  the  first  effort  of  any  subsequent  study. 


Figure  5-1 

84  Flux  Interpo lation  Nodes  of  Hexadic  Three 
Dimensional  Fit 

The  84  polynomials,  which  together  constitute  a  complete 
basis  for  the  hexadic  are  given  in  table  5-1,  with  quantities 
indicated. 


Uith  this  information  the  basis  functions  are  specified. 
The  degrees  of  freedom  F*^^^  and  are  distinct  for 

each  tetrahedron.  and  need  to  be  calculated  only  one 


for  each  triangle  with 


need  to  be  calculated  only  once 


f  — 

(p  =  ^ 


(5-1) 
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Integral  Evaluation 


Chapter  4  results  show  that  the  finite  element  method,  in 
the  case  of  isotropic  scatter  works,  but  one  would  like  to  see 
it  converge  with  less  mesh  refinement,  and  with  a  smaller 
number  of  computations.  A  method  of  exactly  integrating  the 
scattering  terms  is  explained  in  this  chapter  that  should  meet 
this  objective,  as  well  as  restore  the  penalty  function’s 
usefulness  and  guarantee  positive  definite  matrices. 

A.  Hexadic  Interpolation  With  Flux 

A  hexadic  function  in  three  dimensions  requires  84  degrees 
of  freedom  to  be  completely  specified.  If  all  are  flux,  then 
they  can  be  described  in  terms  of  the  twenty  nodal  two 
dimensional  cubic  interpolants  (ten  from  each  triangle), 
independent  of  tetrahedral  geometry.  That  is  to  say,  the  basis 
functions  would  be  constant,  since  they  no  longer  involve 
derivatives  of  natural  coordinates.  There  are  five  distinct 
integrals  to  be  performed,  because  of  the  u  expansion  in  the 
second  scattering  integral.  Basis  functions  can  be  calculated 
separately,  and  stored  in  a  single  matrix  of  dimension  (5,84). 
This  significantly  reduces  calculations,  and  eliminates  the 
requirement  for  the  finite  element  transport  code  to  find  three 
dimensional  interpolants  entirely. 

3.  Interpolation  Nodes  and  Basis  Functions 

The  following  nodes  are  evenly  volume  distributed  and 
should  provide  a  good  hexadic  fit.  Consider  slicing  the 
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element  mesh  D  results  in  nearly  800  tetrahedra  over  which  the 
scattering  integrals  must  be  evaluated.  Each  of  these 
tetrahedra,  treated  in  the  code  as  geometrically  distinct, 
require  separate  interpolation  functions,  and  this  is  the 
probable  cause  of  ca leu  la t iona 1  excesses. 

Summary 

When  scattering  occurs,  the  variational  integral  is 
evaluated  by  integrating  over  space,  angle,  and  scattered  angle 
To  simplify  this  calculation  triangles  are  constrained  in  this 
study  to  columns.  Since  cubic  interpolating  functions  are  used 
for  flux,  the  scattering  integrals  involve  the  product  of  two 
cubics,  or  hexadics.  Mapped  in  three  dimensions,  the  local  and 
nonlocal  triangles  create  tetrahedra. 

Two  methods  were  tried  to  evaluate  these  integrals.  Strict 
numerical  evaluation  with  relations  of  accuracy  up  to  Weddle's 
for  n*6  did  not  obtain  acceptable  accuracy.  Error  with  these 
techniques  was  cumulative,  and  refinement  resulted  in  a  loss  of 
the  global  matrices  positive  definiteness  with  around  15 
triangles  per  mean  free  path,  prior  to  convergence.  Numerical 
evaluation  of  the  scattering  integrals  appears  to  hold  no 
potential  with  the  method's  present  formulation.  Approximating 
the  hexadic  with  another  cubic  gave  better  results,  but  did  not 
guarantee  positive  definiteness.  Solution  accuracy  was 
sufficiently  verified  against  a  Pn  benchmark  over  the  test 
domain.  Convergence  appears  to  occur  with  a  reasonable  but 
larger  amount  of  mesh  refinement  than  in  the  no  scattering  case 
The  number  of  computations  needing  to  be  performed  may  be 
excessive.  Convergence  is  speeded  by  specifying  angle 


FIGURE  4-9 

COMPARISON  OF  Ln  AND  FE  SCALAR  FLUX 
LAMBERTIAN  SOURCE,  WITH  C=5 


finite  element  angular  flux  was  integrated  over  angle  with  the 
assistance  of  IMSL  routine  ICSCCU,  cubic  spline  interpolation 


and 
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A  comparison  of  the  La  and  FE  results  for  mesh  D,  with  a  depth 
of  three  mean  free  paths,  is  displayed  in  table  4-4,  and  graphed 
in  f igure  4-9  . 


X 

0 

.  375 

.75 

1 . 5 

2 .25 

3.0 

La 

.  296 

.191 

.129 

.  064 

C 

.  016 

FE 

.  292 

.  200 

.134 

.061 

.028 

.011 

Table  4-4 

Ln  and  Finite  Element  Comparison  of  Scalar  Fluxes 
Lambertian  Source,  c».5,  Mesh  D,  Fluxes  and  Derivatives 
Specified  as  Boundary  Conditions 

Agreement  between  the  two  codes  is  good.  Differences  are  of 
the  same  order  magnitude  as  the  scattering  error  estimation 
previously  done  for  this  mesh.  At  x»3  ,  the  percentage  difference 
is  large,  but  the  magnitude  of  the  variation  is  small.  The  graph 
of  figure  4-9  displays  the  close  correlation  between  the  two 
separate  calculations  1  results. 

Computationally  the  method  can  be  considered  excessive.  Mesh 
E,  composed  of  46  elements,  requires  over  4  minutes  of  CPU  time 
on  a  Harris  800  computer.  The  correlation  between  Harris  times 
and  the  Vax  times  of  chapter  3  is  unknown,  but  clearly  the 


number  of  calculations  has  greatly  increased.  The  4  column  46 
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considered  as  entirely  due  to  truncation  of  the  polynomial 
series.  Mesh  D  calculations,  with  a  depth  of  three  mean  free 
paths,  and  no  scattering  show  an  average  of  0.127.  difference  at 
u=l  from  analytically  computed  angular  flux.  This  represents  the 
error  from  cubic  approximation  of  flux,  for  the  mesh,  and  degree 
of  refinement  under  scrutiny.  Comparing  this  to  the  scattering 
mesh  D  case  of  three  mean  free  paths  leaves  a  remainder  of 
2.657.,  an  approximate  estimate  of  scattering  integral  evaluation 
error  in  this  case. 

Table  4-2  contains  3  instances  where  less  refined  meshes 
appear  to  give  more  accurate  answers  than  a  denser  mesh.  If 
penalties  were  exact  they  should  indicate,  as  in  the  no 
scattering  case,  that  better  finite  element  fits  can  occur 
without  necessarily  observing  steady  convergence  of  nodal  values 
to  an  "exact"  answer. 

Further  indication  of  the  finite  element  methods  success 
comes  from  investigating  the  lambertian  flux  incident  on  the 
left  boundary  with  c»1.0  .  The  angular  flux  in  this  scattering 
only  medium  of  depth  equal  to  three  mean  free  paths  reflects  the 
hump  predicted  at  x«0  of  figure  4-5.  The  surface  of  angular 
flux,  plotted  in  figure  4-3  shows  that  particles  leak  out  both 
ends,  and  that  angular  flux  is  approaching  isotropy  as  the 
region  is  penetrated. 

Cited  in  Goff's  thesis  (2:67),  were  the  benchmark  case 
results  of  a  transport  code  known  as  Ln .  This  a  program  recently 
developed  as  a  P.H.D.  dissertation  by  LCDR.  Kirk  A.  Mathews 
( AF IT/GNE/ 8 5D ) .  The  output  of  this  code  is  scalar  flux,  so  the 
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CASE  -  integer  reflecting  the  orientation  of  local  and  non 
local  triangles  w.r.t.  each  other 

TIME  -  integer  reflecting  which  half  of  case  2  and  case  4  is 
being  currently  calculated 

PTNODE  -  array  storing  the  global  numbering  of  a  triangles 
finite  element  interpolation  nodes 

COLUMN  -  array  storing  the  column  each  triangle  belongs  to,  the 
top  element  of  that  column,  and  the  number  of  elements 
the  column  posseses 

Variables,  Not  Passed,  by  Subroutine,  Requiring  Definition 
Subroutine  SINFCN 

SGT  -  matrix  of  interpolating  function  constants  for  the 
terahedral  cubic  (2-43) 

M  -  array  storing  the  4x4  partitioned  matrices  of 
(2-42)  and  (2-43) 

Subroutine  SCATA  and  SCATB 

W1 ,  W2  ,  W3  ,W4  ,W5  ,W6  -  dimension  10  vectors  storing  lo'-al  and 
non  local  flux,  and  its  derivatives,  at  locations  that 
are  not  triangular  cubic  interpolation  nodes 
F  -  array  storing  the  twenty  10x10  matrices  used  as 
interpolation  nodes  for  the  three  dimensional  cubic 
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,1,2500 

PROGRAM  FECUBE 

*  FINITE  ELEMENT  SOLUTION  OF  ONE  SPEED  TRANSPORT  EQUATION  IN 

*  SLAB  GEOMETRY,  ISOTROPIC  SCATTER »  CUBIC  APPROXIMATION  OF 

*  FLUX,  CUBIC  APPROXIMATION  OF  HEXADIC  SCATTERING  INTEGRAL 


PARAMETER  <MNQDE=151  ,  MNTRIA=50> 


DOUBLE  PRECISION  CORDND(MNQDE ,2) ,PHI (MNOBE) 


DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
DOUBLE  PRECISION 


AREAS(MNTRIA) ,MA< 10,10) 

SCI <10,33) ,SR1< 18) ,SC2< 10,33) ,SR2<18> 
SC3< 10,33) ,SR3< 18) ,SC4< 10,33) ,SR4<18) 
SC5< 10,33) , SR5< 18) ,SC6< 10,33) ,SR6< 18) 
BC1 (10,20) ,BR1 (10) 

BC2< 10,20) ,BR2< 10) ,D1,D2 
BC3  <10,20), BR3  < 1 0 ) , AS  <  MNODE*  <  MNODE- 1 ) /2 ) 
ML  <  MNTR IA,10,10) 

NLM<MNTRIA, 16, 10,10) ,NLI<MNTRIA,4, 10) 

MG (MNODE, MNODE) 

GT  <  MNTRIA, 10 , 10 ) , LI (MNTRIA, 10 ,4) 

MB( 10,10) ,MS< 10,10) »DRVS( 10,2) 

RANGE , SIGMAT , S I GMAS , PEN (MNTRIA) 
G1,G2,G3,F1,F2,F3,A 
U1,U2,U3»X1 ,X2,X3 
0(5,20) , SGM<5,4,4) 


DOUBLE  PRECISION  E< 4) ,F( 4) ,G< 4) ,H( 5,20) ,SA( 10, 10 ) ,SB( 10 , 10) 
INTEGER  N,TRI, TRIP, NTRIA, CASE, TIME 
INTEGER  PTNODE < MNTRIA ,11), COLUMN (32,2) 

LOGICAL  CHECK1 

COMMON  MG,ML,NLM,NLI ,LI ,GT 


CHECK1  =  .FALSE, 


*  READ  INITIAL  DATA 

CALL  GDATA < NTRIA , N , PTNODE, COLUMN , CORDND , 

C  AREAS , RANGE , SIGMAT , SIGMAS , MA , BC1 , BR1 , BC2 , BR2 , BC3 , BR3 , 

C  SCI , SRI , SC2 , SR2 , SC3 , SR3 , SC4 , SR4 , SC5 , SR5 , SC6 , SR6 , 

C  U,SGM,PHI) 


*  RENUMBER  THE  MESH,  GLOBALLY,  AND  LOCALLY  PER  FIGURE  2-4 
CALL  CHGRID  (PTNODE, CORDND, N, NTRIA) 


*  ZERO  THE  GLOBAL  MATRIX 

DO  69  1=1, N 
DO  68  J»1,N 

MG (I, J) =0.0 

68  CONTINUE 

69  CONTINUE 

*  CALCULATE  PARTICLE  STREAMING, ABSORBING  AND  BOUNDARY  TERMS 

*  ASSEMBLE  INTO  LOCAL  MATRIX  FOR  A  TRIANGLE,  AND  ASSEMBLE 

*  GLOBALLY 

DO  50  TRI=1, NTRIA 

U1=C0RDND< PTNODE <TRI,1), 2) 

A  * 


56  U2=C0RDND< PTNODE< TRI , 4) ,2) 

57  U3=C0RDND< PTNODE < TRI, 7) ,2) 

58  X1=C0RDND<PTN0DE<TRI,1) ,1) 

59  X2=C0RDND<PTN0DE<TRI,4) ,1) 

60  X3=C0RDND< PTNODE (TRI, 7) , 1) 

61  A=AREAS<TRI)*2»0 

62  G1=(U2-U3)/A 

63  G2=<U3-U1)/A 

64  G3=(U1-U2)/A 

65  F1=<X3-X2)/A 

66  F2*(X1-X3)/A 

67  F3=<X2-X1)/A 

68  CALL  INFCN<TRI,G1,G2,G3,F1,F2,F3> 

69  CALL  BNDRY  <U1 ,U2,U3,G1 , G2 ,G3,BC1 , BC2,BC3 , BR1 

70  C  ,BR2,BR3,SIGMAT, MB, DRVS, AREAS, TRI) 

71  CALL  STREAM(SC1,SR1,SC2,SR2,SC3,SR3,SC4,SR4,SC5,SR5, 

72  C  SC6, SR6, MS, U1,U2,U3,G1,G2,G3, AREAS, TRI, 

73  C  DRVS) 

74  CALL  LMATRX  <  MA , MB , MS ,  AREAS ,  S I GMAT ,  TR I  > 

75  CALL  ASEMBL <  PTNCDE ,TRI ) 

76  50  CONTINUE 


78  *  CALCULATE  SCATTERING  CONTRIBUTION  -  FOR  A  TRIANGLE  -  FROM 

79  *  COLUMN  TOP  TO  BOTTOM 

80  DO  150  TRI=1,NTRIA 

81  K=COLUMN< PTNODE (TRI, 11) ,1) 

82  DO  125  TRIP=K,K-1+C0LUMN(PTN0DE(TRI,11),2) 

83  TIME=1 

84  130  CALL  CASEDT ( TRI , TRIP , CORDND , PTNODE ,TIME,E,F,G, 

85  C  V6,CASE,U1,U2,U3,X1 ,X2,X3) 

86  CALL  SINFCN(E,F,G,V,SGM,H) 

87  CALL  SCAT  A ( H , TRI ,  TRIP  , CASE  , TIME , SA ,  CORDND  , PTNODE ) 

88  CALL  SCATB(U1,U2,U3,X1,X2,X3,TRI,TRIP, AREAS, H, 

89  C  , CASE, TIME, SB, CORDND, PTNODE) 

90  CALL  NLMTRX ( TRI , TRIP, SIGMAS, SIGMAT , 

91  C  TIME,U6,SA,SB) 

92  IF  ( CASE •  EQ » 2 , OR  » CASE ♦ EQ ♦ 4 )  THEN 

93  IF  (TIME.EQ.l)  THEN 

94  TIME=TIME+1 

95  GO  TO  130 

96  ENDIF 

97  ENDIF 

98  CALL  SASMBL( PTNODE, TRI, TRIP) 

99  125  CONTINUE 

100  150  CONTINUE 

101 

102  *  PUT  GLOBAL  MATRIX  IN  ITS  QUADRATIC  FORM  - 

103  DO  250  1=1, N 

104  DO  200  J«1,I 

105  MG  <I,J)  =  <MG<I,J) +MG  ( J , I ) ) /2 ► 0 

106  MG<  J, I )*MG< I, J) 

107  200  CONTINUE 

108  250  CONTINUE 

109 

110  *  IF  DESIRED,  DIAGNOSTIC  DATA  CAN  BE  TURNED  ON  IN  ’OUTPUT*  HERE 

111  CALL  OUTPUT (PHI ,N, PTNODE , CORDND, NTR I A, CHECK 1 


, PEN , SIGMAS , RANGE , SIGMAT ) 


112  C 

113  CHECK  1  =*  *  TRUE  ► 

114 

115  *  APPLY  THE  BOUNDARY  CONDITIONS 

116  CALL  BNDCND ( CORDND , PHI , N , NTRIA , RANGE ) 

117 

118  *  PLACE  GLOBAL  MATRIX  IN  BAND  STORAGE  FOR  IMSL 

11?  K»1 

120  DO  350  1*1, N 

121  DO  300  J*1,I 

122  AS( K)=MG( I ,  J) 

123  K=KF1 

124  300  CONTINUE 

125  350  CONTINUE 

126 

127  *  SOLVE  THE  SET  OF  LINEAR  EQUATIONS 

128  CALL  LEQT1P(AS,1,N, PHI, MNODE, IDGT,D1,D2,IER) 

129  PRINT*, 'IER  IS  . ..',IER 

130 

131  *  CALCULATE  PENALTIES,  AND  PRINT  OUT  RESULTS 

132  CALL  PENLTY<PHI,PTNODE, PEN, NTRIA, COLUMN) 

133  CALL  OUTPUT < PHI, N,PTN0DE,C0RDND, NTRIA, CHECK1 

134  C  , PEN, SIGMAS, RANGE, SIGMAT) 

135 

136  END 

137 

138 

139  ************************************************************ 

140 

141  *  GATHER  INITIAL  DATA  -  READS  THREE  DATA  FILES 

142  *  MESH  -  GRID  DATA  (APPENDIX  F) 

143  *  CODATA  -  COEFFICIENTS  OF  LOCAL  MATRICES  (APPENDICES  B,C,D) 

144  *  SDATA  -  CONSTANTS.  FIVE  OF  THE  PARTITIONED  MATRICES  OF  2-42 

145  *  AND  2-43,  AS  WELL  AS  THE  INTEGRALS  OF  BASIS  POLYNOMIALS 

146 

147  SUBROUTINE  GDATA (NTRI A, N,PTNODE, COLUMN, CORDND, 

148  C  AREAS, RANGE, SIGMAT, SIGMAS, MA,BC1,BR1,BC2,BR2,BC3,BR3, 

149  C  SCI , SRI , SC2 , SR2 , SC3 , SR3 , SC4 ,SR4 , SC5 , SR5 , SC6 , SR6 , 

150  C  V,SGM,PHI) 

151 

152  PARAMETER  (MNODE-151  ,  MNTRIA»50> 

153 

154  DOUBLE  PRECISION  CORDND (MNODE, 2) 

155  DOUBLE  PRECISION  AREAS(MNTRIA) ,MA( 10, 10) 

156  DOUBLE  PRECISION  SCI ( 10,33) , SRI ( 18) ,SC2( 10,33) ,SR2( 18) 

157  DOUBLE  PRECISION  SC3( 10,33) ,SR3( 18) ,SC4( 10,33) ,SR4( 18) 

158  DOUBLE  PRECISION  SC5( 10,33) ,SR5( 18) ,SC6( 10,33) ,SR6( 18) 

159  DOUBLE  PRECISION  BC1 ( 10 ,20 ) , BR1 ( 10 ) 

160  DOUBLE  PRECISION  BC2( 10,20) ,BR2( 10) 

161  DOUBLE  PRECISION  BC3( 10,20) ,BR3( 10) 

162  DOUBLE  PRECISION  V(5,20) ,SGM(5,4,4) 

163  DOUBLE  PRECISION  RANGE, SIGMAT, SIGMAS, PHI (MNODE) 

164  DOUBLE  PRECISION  ML(MNTRIA, 10,10) 

165  DOUBLE  PRECISION  NLM(MNTRIA, 16, 10, 10) ,NLI CMNTRIA,4, 10) 

166  DOUBLE  PRECISION  MG ( MNODE , MNODE ) 

167  DOUBLE  PRECISION  GT( MNTRIA, 10 , 10) ,LI(MNTRIA, 10, 4) 
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168  INTEGER  N,NTRIA,TRI ,NB 

169  INTEGER  PTNODE (MNTRI A, 11) ,C0LUMN(32, 2) 

-  170  CHARACTER  TRASH#21 

171  COMMON  MG,ML,NLM,NLI,LI,GT 

172 

173  OPEN  < 15 , FILE* ' MESH ' , STATUS* " OLD ' ) 

174  REWIND  15 

175 

176  READ( 15, ' < A16) ' )  TRASH 

177  READ( 15,'(3(1X,I7>)')  NTRIA,N,NCOL 

178  READ( 13, ' ( IX) ' ) 

179 

180 

181  READ ( 15, ' ( A16) ' )  TRASH 

182  READ(15,'(3(1X,F7.3))')  RANGE, SIGMAT, SIGMAS 

183  RANGE*RANGE*SIGMAT 

184  READ< 13, ' ( IX)  '  > 

185 

186 

187  READ< 15»'(A16)')  TRASH 

188  DO  60  1*1 ,NTRIA 

189  READ(  15, ' ( IX, 17 , 8X,4( IX ,17) ) ' )  TRI , (PTNODE( I , J) , J*1 ,4 ) 

190  60  CONTINUE 

191  READ< 15, '(IX)') 

192 

193  READ < 15,' <A16) ' )  TRASH 

194  DO  70  1*1 ,NCOL 

195  READ< 13,'(3(1X,I7,8X))')  TRI , (COLUMN( I , J) , J*1 ,2) 

A «  196  70  CONTINUE 

197  READ(15,'(1X)') 

198 

199  READ(13,'(A16)')  TRASH 

200  DO  80  1*1, N 

201  READ(13,'(1X,I7,8X,2(2X,F7.3>)'>  NODE, (CORDND( I , J) , J=1 ,2) 

202  CORDND (1,1) *CORDND (1,1) GRANGE 

203  80  CONTINUE 

204  READ< 13, ' ( IX ) ' ) 

205 

206  DO  82  I»1,3*N+NTRIA 

207  PHI < I)*0»0 

208  82  CONTINUE 

209  READ< 15, ' <A16) ' )  TRASH 

210  READ(15»'(I7)'>  NB 

211  DO  83  1*1, NB 

212  READ(15,'(1X,I7,8X,E11.3) ' )  J,PHI(J> 

213  83  CONTINUE 

214 

215  CLOSE  (15) 

216 

217  DO  90  TRI*1 ,NTRIA 

218  U3“C0RDND<PTN0DE<TRI ,3) ,2) 

219  U2*C0RDND  <  PTNODE  <TRI,2) ,2) 

220  X2*C0RDND  <  PTNODE (TRI, 2) ,1) 

221  XI =CORDND( PTNODE ( TRI,1)  ,1) 

222  AREAS ( TRI )*ABS( »3#(U3-U2) *<X2-X1 ) ) 

223  IF  (AREAS(TRI).LT.1.0E-13>  THEN 
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224 

225 

226  90 

227 

228 

229 

230 

231 

232  100 

233 

234 

235 

236 

237 

238 

239  110 

240 

241 

242 

243 

244 

245 

246  120 

247 

248 

249 

250 

251 

252  130 

253 

254 

255 

256 

257 

258 

259 

260  140 

261 
262 

263 

264 

265 

266 

267 

268 

269  150 

270 

271 

272 

273 

274 

275 

276 

277 

278  160 

279 


PRINT**  'AREA  OF  ZERO  IN  ELEMENT''  * TRI 
ENDIF 
CONTINUE 

0PEN(16,FIL£='C0DATA' f STATUS='OLD' ) 

REWIND  16 
DO  100  1=1*10 

READ  (16*'(10(1X*F5*1))')  (MA(I*J)  *J=1*10) 
CONTINUE 


DO  110  1=1*10 

READ  (16*'(10(1X,F5.1))')  <BC1  <  I  ,  J)  *  J»1 , 10) 
READ  (16*/(10(1X*F5*1))/)  < BC1 < I *J> * J=1 1 *20) 
CONTINUE 

READ  (16*'(10(1X»F5*1))')  < BR1 < I ) *  1  =  1  *  10) 

DO  120  1=1,10 

READ  (16*'(10(1X*F5*1>)')  <BC2< I* J) * J=l*10) 

READ  (16*'(10(lX*F5»l))')  (BC2< I *J)*J=11*20) 
CONTINUE 

READ  (16*'(10(1X*F5*1))')  <BR2< I ) *  1=1  *  10) 

DO  130  1=1*10 

READ< 16* ' < 10< 1X*F5. 1 ) ) ' )  ( BC3< I  * J) * J»1 *  10) 
READ< 16* ' ( 10( 1X*F5»1) )  '  )  (BC3( I* J> * J=ll*20) 
CONTINUE 

READ( 16*/(10(1X*F5«1))')  <BR3< I ) *  1=1  *  10 ) 

DO  140  1=1,10 

READ  (16*4200)  <SC1(I* J) * J»l*10) 

READ  (16*4200)  (SCI < I *J> * J=ll *20) 

READ  (16*4200)  (SC1( I  * J) * J=21 *30) 

READ  (16*4100)  (SC1( I* J) * J-31 *33) 

CONTINUE 

READ  (16*4200)  (SRl(I) *1=1* 10) 

READ  (16*4000)  (SRI ( I ) *  1  =  1 1  *  18) 

DO  150  1=1,10 

READ  (16*4200)  ( SC2( I  * J) * J=1 *  10) 

READ  (16*4200)  ( SC2( I  * J ) * J=1 1 *20 ) 

READ  (16*4200)  (SC2( I  * J) * J=21 *30) 

READ  (16*4100)  ( SC2( I  * J) * J=31 *33) 

CONTINUE 

READ  (16*4200)  (SR2( I) *1=1 *10) 

READ  (16*4000)  ( SR2( I ) , 1  =  1 1  *  18 > 

DO  160  1=1*10 

READ  ( 16*420u)  (SC3( I  * J) , J=1 *  10) 

READ  (16,4200)  (SC3( I  * J) * J“ll *20) 

READ  (16*4200)  (SC3( I  * J) * J=21 *30) 

READ  (16*4100)  (SC3( I  * J) * J»31 *33) 

CONTINUE 

READ  (16*4200)  ( SR3< I ) , 1=1 , 10 ) 

A-a 


R£.A,D  ( 16,4000)  (SR3(I) ,1=11,18) 


280 

281 

282 

283 

28A 

285 

286 

287  180 

288 

289 

290 

291 

292 

293 

294 

295 

294  190 

297 

298 

299 

300 

301 

302 

303 

304 

305  200 

306 

307 

308 

309 

310 

311 

312 

313 

314 

315  210 

316 

317 

318 

319  220 

320  230 

321 

322 

323 

324 

325  240 

326  250 

327  260 

328 

329 

330 

331  4000 

332  4100 

333  4200 

334  4300 

335 


DO  180  1=1,10 

READ  (16,4200)  (SC4( I , J) , J=1 , 10) 

READ  (16,4200)  (SC4( I , J) , J=ll ,20) 
READ  (16,4200)  (SC4( I,J) ,J=21,30) 
READ  (16,4100)  ( SC4( I , J) , J-31 ,33) 
CONTINUE 

READ  (16,4200)  (SR4( I ) ,1=1, 10) 

READ  (16,4000)  (SR4( I) , 1=11,18) 

DO  190  1=1,10 

READ  (16,4200)  (SC5( I , J) , J=1 , 10) 
READ  (14,4200)  (SC5( I , J) ,J=11 ,20) 
READ  (16,4200)  (SC5( I , J) , J=21 ,30) 
READ  (16,4100)  ( SC5( I , J) , J=31 ,33) 
CONTINUE 

READ  (16,4200)  ( SR5( I ) , 1=1 , 10) 

READ  (16,4000)  (SRS( I) ,1=11,18) 

DO  200  1=1,10 

READ  (16,4200)  (SC6( I, J) ,J=1, 10) 
READ  (16,4200)  (SC6( I ,J) , J=ll ,20) 
READ  (16,4200)  (SC6( 1, J) ,J=21,30> 
READ  (16,4100)  (SC6(I,J) , J«31,33) 
CONTINUE 

READ  (16,4200)  (SR6( I) , 1=1 , 10) 

READ  (16,4000)  (SR6( I) , 1=11, 18) 

CLOSE  (16) 

OPEN( 17,FILE='SDATA' , STATUS- 'OLD' ) 

REWIND  17 
DO  210  1=1,5 

READ (17,  4200 )  (VK  I, J> , J»i ,10) 

READ( 17,4200 )  ( V( I , J) , J-ll ,20) 
CONTINUE 
DO  230  K=1 ,5 

DO  220  1=1,4 

READ( 17,4300)  (SGM(K , I ,J),J=1,4) 
CONTINUE 
CONTINUE 
DO  260  K=1 , 4 

DO  250  1=1,4 

DO  240  J=l,4 

SGM ( K , I , J ) =SGM (K,I,J)/27»0 
CONTINUE 
CONTINUE 
CONTINUE 
CLOSE  (17) 


FORMAT  (8(1X,F6.1)> 

FORMAT  (3(1X,F6,1)> 

FORMAT  ( 10( IX ,F6 ♦ 1 ) ) 

FORMAT  ( 4( IX ,F6 *  1 ) ) 
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END 


336 

337 

338 

339 

340 

341 

342 

343 

344 

345 

346 

347 

348 

349 

350 

351 

352 

353 

354 

355 

356 

357 

358 

359 

360 

361 

362 

363 

364 

365 

366 

367 

368 

369 

370 

371 

372 

373 

374 

375 

376 

377 

378 

379 

380 

381 

382 

383 
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X  RENUMBERS  THE  GRID  -  SINCE  MESH  IS  NUMBERED  DIFFERENTLY 
X  FOR  EACH  FIT  (LINEAR, QUADRATIC, AND  CUBIC)  ALLOWS  THE  DATA 
X  FILE  *MESH*  TO  REMAIN  SIMPLE,  AND  BE  USED  BY  ALL  THREE 
X  CODES  -  NUMBERING  IS  AS  PER  FIGURE  2-4 

SUBROUTINE  CHGRID  (PTNODE, CORDND, N,NTRIA) 

PARAMETER  (MN0DE“151  ,  MNTRIA*50) 

DOUBLE  PRECISION  CORDND  <  MNODE , 2 ) , D , E 
DOUBLE  PRECISION  ML(MNTRIA, 10 , 10 ) 

DOUBLE  PRECISION  NLM< MNTRIA, 16 , 10 , 10 ) ,NLI < MNTRIA, 4 , 10 ) 

DOUBLE  PRECISION  MG (MNODE, MNODE) 

DOUBLE  PRECISION  GT( MNTRIA, 10 , 10 ), LI ( MNTRIA, 10,4 ) 

INTEGER  N,NTRIA,TRI , A,B,C 
INTEGER  PTNODE( MNTRIA, 11)  . 

COMMON  MG,ML,NLM,NLI ,LI ,GT 


DO  100  I=*l  ,NTRIA 
K“(3XN)+I 

CORDND(K , 1 ) “(  1 » 0/3  »0 )X( CORDND ( PTNODE( 1 , 1 ) , 1 )  +  CORDND ( 
C  PTNODE (I,2),l)  ¥  CORDND (PTNODE( 1,3) ,1) ) 

CORDND (K,2)“(l*  0/3  »  0 )  *  ( CORDND ( PTNODE (1,1), 2)  +  CORDND ( 
C  PTNODE (I, 2), 2)  +  CORDND ( PTNODE (I, 3), 2)) 

100  CONTINUE 


DO  110  TRI“1  ,NTRIA 
A*PTNCDE ( TRI , 1 ) 
B“PTN0DE(TRI,2) 

C“PTNODE( TRI ,3) 

PTNODE ( TRI , 1 1 ) “PTNODE ( TRI , 4 ) 
PTNODE ( TRI , 1 ) “3XA-2 
PTNODE ( TRI , 2 ) *3XA- 1 
PTN0DE(TRI,3)“3XA 
PTNODE (TRI, 4 )»3XB-2 
PTNODE (TRI, 5) “3XB-1 
PTNODE (TRI, 6 )=3XB 
PTNODE( TRI , 7 ) “3XC-2 
PTN0DE(TRI,8)*3XC-1 
PTNODE (TRI, 9 )“3XC 


384  PTNODE (TRI, 10 )=3XN+TRI 


385  110  CONTINUE 

386 

387 

388  DO  120  I“N, 1,-1 

389  D“CORDND( 1,1) 

390  E=CORDND( 1,2) 

391  K“3XI -2 


A— 1 0 


392 

C0RDND(K,1)=D 

393 

CORDND  <  K  *  2 ) *E 

394 

C0RDND(K+1,1)=D 

395 

CORDND (K+l *2)=£ 

396 

CORDND  <  K+2  *  1 )  =*D 

397 

CORDND (KF2 *2 )=£ 

398  120 

CONTINUE 

399 

400 

N»3*N  +  NTRIA 

401 

402 

403 

404 

END 

405 

406 

407 

408 

409  ****%%*t*************************X******X**tt*t**X****************t* 

410 

411  *  FIND  THE  MATRIX  GT,  OF  <2-34) 

412 

413  SUBROUTINE  INFCN<TRI ,G1 ,G2*G3,F1 * F2*F3) 

414 

415  PARAMETER  (MNODE-131  *  MNTRIA-50) 

416  DOUBLE  PRECISION  GT(MNTRIA, 10 * 10 ) *LI < MNTRIA * 10* 4) 

417  DOUBLE  PRECISION  Gl,G2,G3,Fl ,F2,F3,F 

418  DOUBLE  PRECISION  ML (MNTRIA* 10 *10) 

419  DOUBLE  PRECISION  NLM<MNTRIA, 16 » 10* 10) ,NLI < MNTRIA, 4, 10 ) 

iQ*  420  DOUBLE  PRECISION  MG < MNODE * MNODE ) 

421  INTEGER  TRI 

422  COMMON  MG  *  ML  *  NLM  *  NLI *  LI  *  GT 

423 

424  DO  550  1*1*10 

425  DO  300  J»l,10 

426  GT <TRI *  I  * J)*0 .0 

427  300  CONTINUE 

428  530  CONTINUE 

429 

430,  F*G2*F3-G3*F2 

431  GT(TRI*1,1)*1.0 

432  GT  <  TRI , 2 , 1 ) »3*<  63*F 1 -G 1 *F3 ) /F 

433  GT  <  TRI ,2  *2) *F3/F 

434  GT  <TRI ,2,3) »-G3/F 

433  GT(TRI,3*1)»3*<G1*F2-G2*F1)/F 

436  GT (TRI *3»2)*-F2/F 

437  GT (TRI, 3, 3) ®G2/F 

438 

439  F=G3*F1-G1*F3 

440  GT(TRI*4*4)=1»0 

441  GT(TRI,3*4)»3*(G1*F2-G2*F1)/F 

442  GT (TRI *5*3) *F1/F 

443  GT (TRI *3*6)*-Gl/F 

444  GT <  TR I  *  6  *  4 ) =3* <  G2*F3-G3*F2 ) /F 

443  GT (TRI »6*3)*-F3/F 

446  GT ( TRI *6*6) =G3/F 

447 


A-11 


c • 


448  F=G1*F2-G2*F1 

449  GT<TRI,7,7)=1.0 

450  GT ( TR 1 , 8 , 7 ) =3*  <  G2*F3-G3*F2 ) /F 

451  GT(TRI,8,8)=F2/F 

452  GT<TRI,8,9)*-G2/F 

453  GT(TRI,9,7)=3#<G3*F1-G1*F3)/F 

454  GT(TRI,9,8)*-F1/F 

455  GT(TRI,9,9)=G1/F 

456 

457  GT  <  TRI , 10 » 10 ) =27  »0 

458 

459  DO  130  1*1, 7, 3 

460  DO  120  J*I,I+2 

461  GT(TRI,10,I)=GT(TRI, 10 , I >-GT < TRI , J, I ) 

462  GT (TRI, 10,1+1) 3GT < TRI, 10,1+1)  -GT < TRI , J, 1+1 ) 

463  GT<TRI,10,I+2)=GT(TRI,10,I+2)  -GT (TRI , J, 1+2 ) 

464  120  CONTINUE 

465  130  CONTINUE 

466 

467  END 

468 

469 

470  ****************************************************** 

471 

472  *  BOUNDARY  MATRIX  -  ASSEMBLAGE  EXPLAINED  IN  APPENDIX  C 


473 

474 

SUBROUTINE  BNDRY 

(U1 ,U2,U3,G1 ,G2,G3,BC1,BC2,BC3,BR1 

475 

C  , BR2,BR3,SIGMAT,MB,D, AREAS , TR I ) 

476 

477 

PARAMETER  < MNODE- 

>151  ,  MNTRIA*50) 

478 

DOUBLE  PRECISION 

U1 ,U2,U3,G1 ,G2,G3,F 

479 

DOUBLE  PRECISION 

BC1 < 10,20) ,BC2( 10,20) ,BC3( 10,20) 

480 

DOUBLE  PRECISION 

BC( 10,20) 

481 

DOUBLE  PRECISION 

BR1 < 10) ,BR2( 10 ) ,BR3( 10) ,BR( 10) 

482 

DOUBLE  PRECISION 

MB< 10, 10), D( 10, 2) ,SIGMAT ,AREAS(MNTRIA) 

483 

DOUBLE  PRECISION 

ML(MNTRIA,10,10) 

484 

DOUBLE  PRECISION 

NLM  <  MNTRI A ,16,10,10) , NLI (MNTRIA,4, 10) 

485 

DOUBLE  PRECISION 

MG ( MNODE , MNODE ) 

486 

DOUBLE  PRECISION 

GT  <  MNTRI A, 10 ,10) ,LI<  MNTRI A ,10,4) 

487 

INTEGER  TRI 

488 

COMMON  MG,ML,NLM, 

NLI , LI ,GT 

489 

490  * 

ASSEMBLE  THE  DERIVATIVE  MATRIX,  TO  BE  OVERLAYED 

491  * 

NOTE  -  PASSED  AS  DRVS 

492 

D(l,l)=Gl 

493 

D(2,1)=G1 

494 

D(3,1)»G1 

495 

D<4,1)»G2 

496 

D(5,1)»G2 

497 

D<6,t)»G2 

498 

D<7,t)»G3 

499 

D<  8 , 1 ) aG3 

500 

D(9,1)»G3 

501 

D< 10, 1) “G1 

502 

D<1,2)*0.0 

503 

D<2,2)»G2 
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(p* 


304 

D<3,2)*G3- 

505 

D<4,2)»0.0 

506 

D<5,2)*G3 

507 

D<6,2)*G1 

508 

D<7,2)*0.0 

509 

D<  8 ,2) *G1 

310 

D<9,2)=G2 

511 

D< 10,2)*G2 

512 

313  *  MULTIPLY  THE  COEFICIENT  MATRICES  AND  ROWS  BY  APPROPRIATE 
514  *  U  VALUE  -  THEN  SUM 

315  F=SIGMAT*4 . 0*AREAS  <  TRI ) /40320 . 0 


316  DO  100  1-1,10 

517  BR<I)»<U1*BR1<I)+U2*BR2<I)+U3*BR3<I)  >*F 

318  DO  50  J*1 , 20 

31?  BC<I,J)»<U1*BC1<I, J)+U2*BC2<I, J)+U3*BC3<I, J) )  *F 

320  30  CONTINUE 

521  100  CONTINUE 


522 

523  *  OVERLAY  THE  DERIVATIVE  MATRIX  TO  FORM  MB 

524  DO  230  1*1,10 

525  DO  200  J-1,10 

526  MB<I,J)»BC<I,<2*J)-1)*B<I,1)+BC<I,2*J)*D<I,2) 

527  200  CONTINUE 

528  250  CONTINUE 

329 

330  *  AUGMENT  THE  LAST  ROW 

531  DO  300  1*1,10 

532  MB< 10, I)*MB< 10, I )+G3#BR< I ) 

533  300  CONTINUE 


534 

535  *  PLACE  IN  ITS  QUADRATIC  FORM 

536  DO  400  1*1,10 

337  DO  350  J*1,I 

338  MB  <  I ,  J  )  *  <  MB  < I , J ) +MB ( J , I > )/2.0 

539  MB  <  J , I ) *MB  < I , J ) 

340  330  CONTINUE 

541  400  CONTINUE 
342 

543  END 

344 

545 

546  ********************************************************** 


547 

548  *  STREAMING  MATRIX  -  ASSEMBLAGE  EXPLAINED  IN  APPENDIX  D 


34? 

350 

SUBROUTINE  STREAM 

<  SCI , SR 1 , SC2 , SR2 , SC3 , SR3 , SC4 , SR4 , SC5 , SR5 , 

331 

C 

SC6, SR6, MS, U1,U2,U3,G1,G2,G3, AREAS, TRI, 

352 

C 

DRVS) 

353 

334 

PARAMETER  <MNODE= 

131  ,  MNTRIA-50) 

333 

556 

DOUBLE  PRECISION 

AREAS<MNTRIA) 

537 

DOUBLE  PRECISION 

SCI < 10,33) , SRI ( 18) ,SC2( 10,33) ,SR2<18> 

338 

DOUBLE  PRECISION 

SC3< 10,33) ,SR3( 18) ,SC4< 10,33) ,SR4< 18) 

339 

DOUBLE  PRECISION 

SC3< 10,33) ,SR5<19) ,SC6( 10,33) ,SR6<18) 

Afl  3 


560 

DOUBLE 

PRECISION 

GG<3) , DSC  10,33) , SC (10,33) ,A,B,C,D,E,F,G 

561 

DOUBLE 

PRECISION 

SR< 18) , DR< 18) 

562 

DOUBLE 

PRECISION 

MS (10, 10) , DRVS  <10,2) 

563 

DOUBLE 

PRECISION 

G1  ,G2 ,G3 ,U1 ,U2 ,U3 

564 

DOUBLE 

PRECISION 

ML(MNTRIA, 10,10) 

565 

DOUBLE 

PRECISION 

NLM<MNTRIA, 16,10,10) ,NLI(MNTRIA,4,10) 

566 

DOUBLE 

PRECISION 

MG  <  MNODE , MNODE  > 

567 

DOUBLE 

PRECISION 

GT (MNTRIA, 10 ,10) ,LI(MNTRIA,10,4) 

568 

INTEGER  TRI 

569  COMMON  MGpML, NLM,NLI ,LI ,GT 

570 

571  *  ASSEMBLE  THE  MATRIX  OF  DERIVATIVES 

572  GG<1)»G1 

573  GG(2)*G2 

574  GG  <  3  )  *G3 

575 


576  *  FILL  IN  COLUMNS  1,2,3,12,13,14,23,24,25  OF  DS 

577  DO  110  1*1, 7, 3 

578  L*l+< 1-1 >/3 


579 

580 

581 

582 

583 

584  100 


K»l  +  <  <I-l)*ll/3> 

DO  100  J*l,10 

DS  <  J , K  >  *DRVS ( J , 1 ) #GG ( L ) 

DS  <  J  ,  K+ 1 )  *DR  VS  <  J  ,  2  )  XGG  (  L  ) 
DS<  J,K+2) *0 » 0 
CONTINUE 


585  110  CONTINUE 

586 


587 
%  588 

589 

590 


DS<10,3)*G1XG3 
DS(10,14)*G2*G3 
DS< 10,25)*G3XG3 


591  *  FILL  IN  REMAINING  COLUMNS 

592  DO  120  J»l,10 

593  DS< J,4)»DRVS< J,1)*G1 

594  DS<  Jp5) *DRVS<  J,2) *G1 

595  DS< J,6)»DRVS<J,1)*G2 


596 

597 

598 

599 

600 
601 
602 


DS  <  J , 7 ) *DRVS ( J , 2 ) XG2 
DS  ( J  ,  8 ) *DS  <  J , 4 ) 

DS  <  J  ,  9 )  *wS  (  J  ,  5  ) 

DS< J, 10)=DRVS< J,1)*G3 
DS  <  J , 1 1 ) = DRVS <  J , 2 ) XG3 
DS( J, 15)=DS< J,6) 

DS ( J  ,  16 ) *DS < J,7) 


603 

DS<  J, 17)*DS<  J, 10 ) 

604 

DS ( J , 18 ) *DS <  J , 1 1 ) 

605 

DS< J,19)*DS< J,6) 

606 

DS  <  J , 20 ) *DS ( J , 7 ) 

607 

DS<  J,21 )*DS<  J,4) 

608 

DS<  J ,22 ) *DS < J , 5 ) 

609 

DS< Jr26)*DS( J,10) 

610 

DS  ( J  ,  27 )  *DS  <  J ,  1 1 ) 

611 

DS ( J , 28 ) *DS  <  J , 4 ) 

612 

DS  ( J  ,  29 )  *DS  <  J  ,  5 ) 

613 

DS  (  J , 30 ) *  DS  <  J , 1 0 ) 

614  DS<J,31)*DS< J,ll) 

615  DS< Jp32)=DS< J,6) 
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616  DS< J*33)=DS< J*7) 

617  120  CONTINUE 

618 

619  DS<10,6)=G1*G3 

620  DS< 10*7) =0*0 

621  DS< 10  *  10 )=G1*G3 

622  DS( 10  *  11 ) =0  *0 

623  DS (10*17) =G2#G3 

624  DS( 10* 18)=0  *0 

625  DS<10,21)=G2*G3 

626  DS< 10*22)=0*0 

627  DS<10*2B)=G3*G3 

628  DS< 10*29)=0  *0 

629  DS< 10*32)=G3*G3 

630  DS<10*33)=0»0 

631 

632  0R(1)=G2*G2 

633  DR(2)=G2*G3 

634  DR(3)=DR(2) 

635  DR<4)=G3*G3 

636  DR<5)»G1*G3 

637  DR<6) =DR(  4) 

638  DR<7)«G1*G1 

639  DR<8)=DR<5) 

640  DR<9)=DR<7) 

641  DR<10)=G1*G2 

642  0R< 11 )=DR< 10 ) 

643  DR< 12) =DR< 1 ) 

644  DR<13)=DR(7) 

645  DR< 14)=DR< 10) 

646  DR<15)»DR<5) 

647  DR< 16)=DR< 1 ) 

648  DR<17)=DR<2) 

649  DR( 18)=DR<  4) 

650 

651  *  MULTIPLY  THE  COEFICIENT  MATRICES  AND  ROUS  BY  THE 

652  *  APPROPRIATE  U'S  -  THEN  SUM 

653  A»U1*U1 

654  B-U2*U2 

655  C=U3*U3 

656  D=Ul*U2*2.0 

657  E=U2*U3#2.0 

658  F=Ui*U3*2,0 

659  G=2 . OfcAREAS  <  TRI ) /40320  » 0 

660  DO  140  1=1*10 

661  DO  130  J=1 ,33 

662  SC(I* J)=<SC1 < I* J)*A  +  SC2< I*J)tB  +  SC3(I,J)*C  + 

663  C  SC4<I,J)*D  ¥  SC5< I  * J)*E  +  SC6 < I  * J ) *F ) *G 

664  130  CONTINUE 

665  140  CONTINUE 

666  DO  150  1=1*18 

667  SR< I ) *<SR1 < I ) *A  ¥  SR2<I)*B  +  SR3<I)*C  + 

668  C  SR4<I)*D  +  SR5<I)*E  +  SR6<I)*F)*G 

669  150  CONTINUE 

670 

671  *  COMPUTE  COLUMNS  1*4,7  OF  STREAMING  MATRIX 

^15 


c 


672 

673 

674 

675 

676 

677 

678 

679 

680 
681 
682 

683 

684 

685 

686 

687 

688 

689 

690 

691 

692 

693 

694 

695 

696 

697 

698 

699 

700 

701 

702 

703 

704 

705 

706 

707 

708 

709 

710 

711 

712 

713 

714 

715 

716 

717 

718 

719 

720 

721 

722 

723 

724 

725 

726 

727 


160 

170 


DO  170  1=1,10 

DO  160  J=1 ,7,3 

K=1  +  (  ( J-l)*U/3> 

MS<I,J)=SC<I,K)*DS<I,K)  +  SC<I,K+1)*DS<I  K+l) 
+  SC<I,K+2)*DS<I,K+2) 

CONTINUE 

CONTINUE 


*  COLUMNS  2,5,  AND  8 

DO  190  1=1,10 

DO  180  J=2 ,8,3 

K=4+<  <  J-2)  *11/3) 

MS  < I , J ) =SC  < I , K ) #DS ( I , K )  +  SC ( I ,K+1 ) #DS< I , K+l ) 

C  +  SC<I,K+2)*DS<I,K+2)  +  SC< I ,K+3) *DS< I ,K+3) 

K=K+4 

MS<I,J+1)=SC<I,K)*DS<I,K)  +  SC( I ,K+1 )#DS( I ,K+1 ) 

C  +  SC( I , K+2)#DS< I ,K+2)  +  SC < I ,K+3> *DS ( I ,K+3) 

180  CONTINUE 

190  CONTINUE 

*  AUGMENT  THE  LAST  ROW 

DO  200  1=2, 8, 3 

K=1  +  4*<I-2)/3 

MS( 10,1) *MS< 10,1) +SR<  K ) *DR<  K ) +SR (K+l >  #DR(K+1 ) 
MS<10,I+l)=MS(10,I+l)+SR<K+2>*DR<K+2)+SR<K+3)*DR<K+3) 
200  CONTINUE 


MS< 10, 10)=0.0 
DO  210  1=13,18 

MS< 10, 10 )=MS< 10 , 10)+SR  < I )#DR< I ) 
210  CONTINUE 


*  FORM  COLUMN  JO  WITH  SYMMETRY 
DO  220  1=1,9 

MS< I , 10 )=MS< 10,1) 

220  CONTINUE 

END 


**********************;»**********:M**************************** 

*  LOCAL  MATRIX  -  MULTIPLY  ABSORBING, BOUNDARY  AND  ~TREAMING 
t  BY  APPROPRIATE  CONSTANTS  -  AND  SUM 

*  PRE  AND  POST  MULTIPLY  BY  GT  TO  COMPLETELY  PRT^ARE  FOR 

*  GLOBAL  ASSEMBLAGE 


SUBROUTINE  LMATRX < MA , MB , MS , AREAS , S I GMAT , TR I ) 

PARAMETER  <MN0DE=151  ,  MNTRIA=50 ) 

DOUBLE  PRECISION  AREAS (MNTRIA) ,MA< 10, 10) 

DOUBLE  PRECISION  M9( 10 , 10 ) ,MS( 10 , 10 ) ,MT < 10, 10 ) 

DOUBLE  PRECISION  SIGMAT,F 
DOUBLE  PRECISION  ML<MNTRIA, 10, 10) 

DOUBLE  PRECISION  NLM< MNTRIA, 16,10,10) ,NLI( MNTRIA ,4,10) 
DOUBLE  PRECISION  MG < MNODE , MNODE ) 

DOUBLE  PRECISION  GT < MNTRIA , 10 , 10 ) , LI ( MNTRIA , 10 , 4 ) 
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UNCLASSIFIED 


NL 


MICROCOPY  RESOLUTION  TEST  CHART 

NAIrnhHl  buht-AU  Of  STANDARD^  TQM  A 


6* 


728 

729 

730 

731 

732 

733 


INTEGER  TRI 

COMMON  MG,ML,NLM,NLI,LI,GT 
F=SIGMAT*SIGMAT*2 ♦ 0#AREAS  <  TRI ) /40320 . 0 


734 

735 

736 

737 

738 

739 

740 

741 

742 

743 

744 

745 

746 

747 

748 

749 

750 

751 

752 

753 

754 

755 

756 

757 

758 

759 

760 

761 


DO  830  1=1,10 

DO  820  J=1 ,10 

ML  <  TRI , J , I ) =MB  <  J  »  I )  +  MA<J,I)*F  +  MS<J,I> 

820  CONTINUE 

830  CONTINUE 

DO  860  1=1,10 

DO  850  J=l,10 
MT<I»J)=0,0 
DO  840  K=1 , 10 

MT  < I , J  >  =MT  < I , J  >  +  GT<TRI,K,I)*ML<TRI,K,J) 

840  CONTINUE 

850  CONTINUE 

860  CONTINUE 

DO  890  1=1,10 

DO  880  J=1 , 10 

ML < TRI ,1,J)=0»0 
DO  870  K=1 , 10 

ML < TRI , I, J)=ML(TRI , I, J)  +  MT< I ,K)*GT(TRI ,K, J) 
870  CONTINUE 

880  CONTINUE 

890  CONTINUE 

END 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


762  *  ASSEMBLE  LOCAL  TERMS  GLOBALLY 

763 


764 

765 

766 

767 

768 

769 

770 

771 

772 
77  3 

774 

775 

776 

777  900 

778 

779 

780 

781 

782  910 

783  920 


SUBROUTINE  ASEMBL<PTNODE,TRI ) 

PARAMETER  <MN0DE=151  ,  MNTRIA-50) 

• 

DOUBLE  PRECISION  ML<MNTRIA, 10 , 10 ) 

DOUBLE  PRECISION  NLM< MNTRIA, 16 , 10, 10) ,NLI < MNTRIA,4, 10) 
DOUBLE  PRECISION  MG <  MNODE , MNODE ) 

DOUBLE  PRECISION  GT< MNTRIA, 10, 10), LI< MNTRIA, 10, 4) 
INTEGER  PTN0DE<MNTRIA,11) ,TRI,R<10) 

COMMON  MG,ML,NLM,NLI ,LI ,GT 

DO  900  1=1,10 

R ( I ) =PTNODE  <  TR I , I ) 

CONTINUE 

DO  920  1=1,10 

DO  910  J*l,10 

MG<R<I) ,R< J) )=MG<RCI) ,R< J) )  +  ML(TRI,I,J) 
CONTINUE 

CONTINUE  A-17 


^  ft  ft  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  j^C  ft  ^  ^  ^  3^C  ^C  ^C  ^C  ^C  ^C  ^fC  ft  ^  ^  ^  ^  ^  ^  ^  ||^ 

*  INSURE  FLUXES  (AND  IN  THIS  CASE  U-CURRENTS)  ARE  AS  SPECIFIED 

*  ON  BOUNDARIES  -  IF  FLUXES  ONLY  ARE  TO  BE  SPECIFIED  DELETE 

*  THE  1+2  TERMS  MODIFICATION 

SUBROUTINE  BNDCND < CORDND , PHI , N , NTRI A , RANGE ) 

PARAMETER  (MN0DE=151  ,  MNTRIA=50> 

DOUBLE  PRECISION  CORDND <  MNODE , 2  > , PHI ( MNODE > 

DOUBLE  PRECISION  ML(MNTRIAf10,10> 

DOUBLE  PRECISION  NLM< MNTRIA, 16 , 10 , 10) ,NLI (MNTRIA, 4 , 10) 
DOUBLE  PRECISION  MG < MNODE  , MNODE  > 

DOUBLE  PRECISION  GT(MNTRIA, 10 , 10 ), LI (MNTRIA, 10,4 ) 

DOUBLE  PRECISION  RANGE 

INTEGER  N, NTRI A 

COMMON  MG,ML,NLM,NLI,Lr,GT 

DO  120  1=1 ,N— NTRIA,3 

IF  <CORDND<Ifl) .EQ*O.O.AND»CORDND(I,2) .GE»0»0>  THEN 
MG<I,I)=MG<I,I)*l*0E+20 
PHI(I)=PHI(I)*MG(I,I) 

MG < 1+2 , 1+2 ) »MG < 1+2 f 1+2 > *1 . 0E+20 
PHI < 1+2) =PHI < 1+2) *MG  < 1+2 , 1+2  > 

ELSE 

IF  <C0RDND<I,1> .EQ.RANGE,AND,C0RDND<If2) ,LE.0,0>  THEN 
MG(I,I)=MG<I,I)*l,0E+20 
PHI<I)=PHI<I)#MG<I,I> 

MG< 1+2  f 1+2 ) =MG < 1+2 , 1 +2  >  *1 , 0E+20 
PHI < I+2)=PHI < I+2)*MG< 1+2, 1+2) 

ENDIF 

ENDIF 

120  CONTINUE 


***:M'************************)M************3Mc******************* 
*  CALCULATE  VALUE  OF  VARIATIONAL  INTEGRAL  OVER  AN  ELEMENT 


SUBROUTINE  PENLTY ( PHI , PTNODE f PEN , NTRI A  f COLUMN ) 

PARAMETER  <MN0DE=151  ,  MNTRIA=50 ) 

DOUBLE  PRECISION  ML < MNTRIAfIO f 10 ) fPHI (MNODE) fPEN< MNTRIA) 
DOUBLE  PRECISION  P<10),L<10) 

DOUBLE  PRECISION  S<10)fF 

DOUBLE  PRECISION  NLM< MNTRIA, 16 f 10 f 10 ) ,NLI ( MNTRIA, 4 , 10 ) 
DOUBLE  PRECISION  MG (MNODE, MNODE) 

DOUBLE  PRECISION  GT< MNTRIA, 10, 10 ), LI (MNTRIA, 10 , 4 ) 

INTEGER  TRI , NTRIA , PTNODE ( MNTRIA , 1 1 ) 


INTEGER  TRIP, COLUMN <32,2) 

COMMON  MG,ML,NLM,NLI,LI,GT 

DO  100  TRI=1, NTRIA 

*  LOCAL  MATRIX  CONTRIBUTION 

DO  20  1=1,10 

P  < I ) =PHI (  PTNODE ( TRI , I ) ) 

20  CONTINUE 

DO  40  1=1,10 
L< I )=0  *0 
DO  30  J=1 , 10 

L< I  )=L( I )  +  ML < TRI »I,J)#P(J) 

30  CONTINUE 

40  CONTINUE 

PEN<  TRI ) =0  »0 
DO  50  1=1,10 

P£N<TRI)=P£N<TRI )  +  L<I)*P(I) 

50  CONTINUE 

*  SUM  OF  NON  LOCAL  MATRICES  CONTRIBUTIONS 

K=COLUMN  <  PTNODE  <  TR 1 , 1 1 ) , 1 ) 

DO  70  TRIP=K,K-1+C0LUMN( PTNODE (  TRI , 11) ,2) 
DO  55  1*1,10 

S< I)*PHI<PTNODE<TRIP,I>  > 

55  CONTINUE 

DO  65  1*1,10 
L(I)=0.0 
DO  60  J»l,10 

L<I)=L<I)  +  NLM<TRI,TRIP,I,J>*S<J> 
60  CONTINUE 

65  CONTINUE 

DO  68  1=1,10 

PEN  ( TRI ) =PEN (TRI)+L(I)*P(I) 

68  CONTINUE 

70  CONTINUE 

PEN  <  TR I ) * . 5*PEN  <  TRI ) 

100  CONTINUE 


*##*«*****»**»#***»*****#»#*******«*******************#****#**** 

*  PRINT  OUTPUT  -  COMPARE  FE  SOLUTION  WITH  Pn  OR  ANALYTICAL 

*  IF  NO  SCATTER 

SUBROUTINE  OUTPUT < PHI , N , PTNODE , CORDND , NTRIA , 

C  CHECK 1 , PEN , S I GMAS , RANGE , S I GM AT ) 


PARAMETER  (MNODE-151  ,  MNTRIA=50> 

DOUBLE  PRECISION  PHI<MNODE) , CORDND (MNODE, 2) 

DOUBLE  PRECISION  PEN (MNTRIA) ,PENTOT 
DOUBLE  PRECISION  ML<MNTRIA, 10, 10) 

DOUBLE  PRECISION  NLM<MNTRIA,16,10,10) ,NLI(MNTRIA,4,10) 
DOUBLE  PRECISION  MG (MNODE, MNODE) 

DOUBLE  PRECISION  GT (MNTRIA, 10, 10) ,LI (MNTRIA, 10,4 ) 
DOUBLE  PRECISION  RANGE , SIGMAT , TPEN 
INTEGER  PTNODE (MNTRIA, 11) , TRI, N, NTRIA 
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LOGICAL  CHECK 1 

COMMON  MG,ML,NLM,NLI,LI,GT 


896 

897 

898 

-  899  IF  < CHECK 1)  THEN 

900 

901  PRINT*, 'NTRIA  N  SIGMAS' 

902  WRITE  (*,4999)  NTRIA,N,SIGMAS 

903  4999  FORMAT  <3X, 13, 5X , 13, 5X ,F6 .3) 

904  PRINT*, 'RANGE  IS*  »  » . ' , RANGE 

905 

906  PRINT*, 'NODAL  VALUES  OF  THE  FLUX' 

907  J=N-NTRIA 

908  DO  100  1*1, J, 6 

909  K=l+< 1-1 )/3 

910  WRITE<*,6010)  K,PHI< I) ,K+1 ,PHI < 1+3) 

911  6010  FORMAT  <2(2X, 13, 3X ,F9  *4) ) 

912  100  CONTINUE 

913 

914  PRINT#, 'ELEMENT  PENALTY  VALUES' 

915  PENTQT=0.0 

916  TPEN=0.0 

917  DO  110  1*1, NTRIA, 2 

918  WRITE<*,6221)  I,PEN<I) ,I+1,PEN<I+1) 

919  PENTOT-PENTOT  +  ABS<PEN<I>)  +  ABS<PEN< 1+1 > > 

920  TPEN-TPEN  +  PEN<I)  +  PEN(I+1) 

921  110  CONTINUE 

922  PRINT*, 'TOTAL  PENALTY  AND  SUM  OF  ABS< PENALTY)  ARE 

923  WRITE <*,6222)  TPEN,PENTOT 

a #  924  6221  FORMAT  <2<2X,I3,5X,E11*5) ) 

925  6222  FORMAT  <2< 10X,E11.5> ) 

926 

927  *  COMPARE  FE  SOLUTION  WITH  APPROPRIATE  BENCHMARK 

928  IF  < SIGMAS, EQ. 0.0)  THEN 

929  CALL  ANALY<PHI,CORDND,SIGMAT,N, NTRIA, RANGE) 

930  ELSE 

931  CALL  PN<PHI,CORDND, SIGMAS, N, NTRIA, RANGE) 

932  END  IF 

933 

934  ELSE 

935  *  IF  DESIRED  TURN  ON  DIAGNOSTIC  OUTPUT  HERE 

936  GO  TO  301 

937  PRINT*,  'MESH  DEFINITION' 

938  PRINT*,  '  TRIANGLE  GLOBAL  NODES' 

939  DO  50  TRI=1, NTRIA 

940  WRITE<*,6050)  TRI , <PTNODE< TRI , I > , 1=1 ,7,3) 

941  6050  FORMAT  <4X,I3,8X,3<I3»3X)) 

942  50  CONTINUE 

943 

944  PRINT# 

945  PRINT*,'  NODE  COORDINATES  <X,U)' 

946  DO  60  I»1,N-NTRIA,3 

947  K=< 1+2) /3 

948  WRITE< *,6060)  K , (CORDND< I,J),J=1,2) 

949  6060  FORMAT <4X,I3»7X»2<F7.3,3X) ) 

950  60  CONTINUE 

951 


A-  20 


932  PRINT* 

933  PRINT*, 'GLOBAL  MATRIX' 

934  PRINT*. 

955  DO  300  1=1, N 

956  WRITE<*,6210)  <MG< I , J) , J=1 ,N) 

937  6210  FORMAT <lX,*t*,16<lX,F6*3) ) 

938  URITE(*,6220) 

959  6220  FORMAT < IX, ‘J •) 

960  300  CONTINUE 

961  301  END IF 

962 

963  350  END 

964 

965  *********************************************************** 

966 

967  *  ASSEMBLE  SCATTERING  MATRICES  <NON  LOCAL)  -  A  SEPARATE 

968  *  SUBROUTINE  IS  USED  BECAUSE  DIMENSIONS  OF  NLM  ARE  DIFFERENT 

969  *  THAN  ML 

970 

971  SUBROUTINE  SASMBL<PTNOBE,TRI ,TRIP) 

972 

973  PARAMETER  <MN0DE>151  ,  MNTRIA-50) 

974  DOUBLE  PRECISION  ML(MNTRIA, 10, 10) 

975  DOUBLE  PRECISION  NLM<MNTRIA, 16, 10, 10) ,NLI(MNTRIA,4, 10) 

976  DOUBLE  PRECISION  MG ( MNODE , MNODE ) 

977  DOUBLE  PRECISION  GT(MNTRIA,10,10) ,LI(MNTRIA,10,4) 

978  INTEGER  PTN0DE<MNTRIA,11) ,TRI,RC10)  ,TRIP 

979  INTEGER  L<10) 

980  COMMON  MG, ML, NLM, NLI, LI ,6T 

981 

982  DO  900  1-1,10 

983  R< I)«PTNODE<TRI , I) 

984  L  < I ) -PTNODE  <  TRIP , I ) 

983  900  CONTINUE 

986 

987  DO;  920  1=1,10 

988  DO  910  J»l,10 

989  MG(R( I ) ,L<  J) )=MG<R< I ) ,L<  J) )  +  NLM<TRI,TRIP,I, J) 

990  910  CONTINUE 

991  920  CONTINUE 

992  END 

993 

994 
993 

996  ************************************************************ 

997 

998  *  COMPARE  FE  SOLUTION  TO  ANALYTICAL  IN  THE  CASE  OF  NO  SCATTER 

999 

1000  SUBROUTINE  ANALY<PHI ,CQRDND,SIGMAT ,N,NTRIA,RANGE) 

1001 

1002  PARAMETER  <MN0DE*131  ,  MNTRIA-50) 

1003  DOUBLE  PRECISION  C0RDND<MN0DE,2) ,PHI(MN0DE) , A, RANGE 

1004  DOUBLE  PRECISION  PERC,TPERC 

1005  DOUBLE  PRECISION  ML<MNTRIA, 10, 10) 

1006  DOUBLE  PRECISION  NLM< MNTRIA, 16, 10 , 10) ,NLI <MNTRIA,4, 10 ) 

1007  DOUBLE  PRECISION  MG (MNODE, MNODE) 


1008 

1009 

1010 
1011 
1012 

1013 

1014 

1015 

1016 

1017 

1018 

1019 

1020 
1021 
1022 

1023 

1024 

1025 

1026 

1027 

1028 

1029 

1030 

1031 

1032 

1033 

1034 

1035 

1036 

1037 

1038 

1039 

1040 

1041 

1042 

1043 

1044 

1045 

1046 

1047 

1048 

1049 

1050 

1051 

1052 

1053 

1054 

1055 

1056 

1057 

1058 

1059 

1060 
1061 
1062 
1063 


DOUBLE  PRECISION  GT< MNTRIA, 10 , 10 > *LI < MNTRIA, 10,4 ) 

INTEGER  N,NTRIA 

COMMON  MG,ML,NLM,NLI,LI,GT 

TPERC=0.0 

K=0 

PRINT#*'  COORDINATES  CURRENTS  FIN  ELEM' 

PRINT#* '  X  U  X  U  FLUX  FLUX  7.  BIFF' 

DO  100  1*1  *  N-NTRI A ,  3 

IF  <CORDND( 1*2). GT .0.0)  THEN 

A-CORDND  <1*2  >*EXP  < -S I GMAT/CORDND < I  *  2 ) *CORDND < I  *  1 ) ) 
PERC»100*ABS  <  PHI < I > -A) /A 

URITE<#*5002)  CORDND< I *1 ) *CORDND< 1*2) *PHI < 1+1 ) *PHI < 1+2)  . 

C  *  PHI ( I ) *  A  *  PERC 

TPERC*TPERC+PERC 

IF  <C0RDND<I*1) .NE. 0.0. AND. C0RDND<I,1) .NE. RANGE)  THEN 
K-K+l 
END  IF 
ENDIF 

100  CONTINUE 

PRINT#, 'AVERAGE  X  DIFFERENCE  IS  ..',TPERC/K 
D»NTRIA*. 5/RANGE 

PRINT*, 'FOR  AN  AVG.  OF' ,D, 'TRIANGLES  PER  MFP  FOR  U>0  ' 

5002  F0RMAT<6<2X*F6.3) *2X*F6.2) 

END 

X*#*********##**###***#****#***#**#******#*******#**#***#******#* 

*  CALCULATE  THE  NON  LOCAL  MATRIX  FOR  TRIANGLE  TRI 

*  INTO  TRIANGLE  TRIP 

SUBROUTINE  NLMTRX (TRI, TRIP* SIGMAS, SIGMAT* 

C  TIME*V6*SA*SB> 

PARAMETER  (MNODE-151  ,  MNTRIA-50) 

DOUBLE  PRECISION  ML<MNTRIA* 10* 10 ) 

DOUBLE  PRECISION  NLM<MNTRIA* 16* 10*10) »NLI<MNTRIA*4* 10) 

DOUBLE  PRECISION  MG < MNODE * MNODE ) 

DOUBLE  PRECISION  GT< MNTRIA, 10, 10 ) *LI (MNTRIA, 10*4 ) 

DOUBLE  PRECISION  SIGMAS*SIGMAT*SA( 10*10) *SB< 10*10) 

DOUBLE  PRECISION  V6 
INTEGER  TRI *TRIP*TIME 
COMMON  MG*ML*NLM*NLI*LI,GT 

*  CALCULATE  CONSTANTS 

A* » 5*S IGMAS*S I GMAS  -  SIGMAS#SIGMAT 

B—SIGMAS 

F*V6#A/720.0 

G*V6#B/5040 . 0 

*  ZERO  THE  NON  LOCAL  MATRIX 

IF  (TIME.EQ.l)  THEN 
DO  50  1*1,3 

DO  40  J*1 *3 

NLM  <TRI,TRIP*I,J)*0.0 
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1064  40  CONTINUE 

1065  50  CONTINUE 

1066  ENDIF 

1067 

1068  *  CALCULATE  THE  NON  LOCAL  MATRIX 

1069  DO  100  1*1,10 

1070  DO  60  J*l,10 

1071-  NLM<TRI,TRIP,I, J)=NLM<TRI,TRIP,I, J)+<F*SA< I, J) 

1072  C  +G*SB<I,J>> 

1073  60  CONTINUE 

1074  100  CONTINUE 

1075 

1076  END 

1077  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

1078  *  FIND  PHI  OF  <L1,L2,L3)  FOR  THE  TRIANGLE  IN  QUESTION 

1079 

1080  SUBROUTINE  PHII<TRI ,L1 ,L2,L3,D) 

1081 

1082  PARAMETER  <MN0DE=151  ,  MNTRIA=50) 

1083  DOUBLE  PRECISION  LI ,L2,L3, D< 10) ,U< 10) 

1084  DOUBLE  PRECISION  ML(MNTRIA, 10, 10) 

1085  DOUBLE  PRECISION  NLM<MNTRIA, 16, 10, 10) ,NLI<MNTRIA,4,10) 

1086  DOUBLE  PRECISION  MG ( MNODE , MNODE ) 

1087  DOUBLE  PRECISION  GT<MNTRIA , 10 , 10) ,LI CMNTRIA, 10 ,4) 

1088  INTEGER  TRI 

1089  COMMON  MG,ML,NLM,NLI,LI,GT 

1090 

1091  W<1)»L1**3 

1092  U<2)»L2*L1**2 

1093  W<3)-L3«L1**2 

1094  U<4)»L2**3 

1095  <4<5)»L3*L2**2 

1096  U<6)»L1*L2**2 

1097  U<7)=L3**3 

1098  U<8)»L1*L3**2 

1099  U<9)»L2*L3**2 

1100  U<10)»L1*L2*L3 

1101  DO  30  1*1,10 

1102  D<I)»0*0 

1103  DO  20  J=l,10 

1104  D(I)=D<I)+W< J)*GT<TRI,J,I) 

1105  20  CONTINUE 

1106  30  CONTINUE 

1107 

1108  END 

1109  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

1110 

1111  *  FIND  D<PHI)/DX  FOR  THE  TRIANGLE  UNDER  SCRUTINY 

1112 

1113  SUBROUTINE  PHIXCTRI ,L1,L2,L3,G1 ,G2,G3,DX) 

1114 

1115  PARAMETER  <MN0DE=151  ,  MNTRIA=SO) 

1116  DOUBLE  PRECISION  LI ,L2,L3,G1 ,G2,G3,W< 10) ,DX< 10) 

1117  DOUBLE  PRECISION  ML<MNTRIA, 10, 10) 

1118  DOUBLE  PRECISION  NLM(MNTRIA, 16, 10, 10) ,NLI <MNTRIA,4, 10) 

1119  DOUBLE  PRECISION  MG < MNODE, MNODE) 
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1120  DOUBLE  PRECISION  GT< MNTRIA , 10 , 10> , LI ( MNTRIA, 10 , 4 ) 

1121  INTEGER  TRI 

1122  COMMON  MG,ML,NLM,NLI,LI,GT 

1123 

1124 

1125 

1126  *  FIND  D<PHI)/DX 


1127  I4<1)=3.0*G1*L1**2 

1128  14  ( 2 )  =G2*L1 **2+L2*Gl*2  »  0*L1 


1129  U<3)=G3*L1**2-H-3*G1*2.0*L1 


1130 

1131 

1132 

1133 

1134 

1135 

1136 

1137 

1138 

1139 

1140 

1141  40 

1142  50 

1143 

1144 


U<4)=3.0*G2*L2**2 
U  <  5 )  =G3*L2**2+L3*G2*2 . 0  *L2 
U <  6 ) =G1*L2**2+Li *G2*2 . 0*L2 
14<7)=3.0#G3*L3**2 
U <  8 ) =G1 *L3**2+L1 *G3*2 .  0*L3 
W<  9 ) =G2*L3**2+L2*G3*2  ,0*L3 
U < 10 ) =G1*L2*L3+G2*L 1 *L3+G3*L1 *L2 
DO  50  1=1,10 
DX< I )=0.0 
DO  40  J=l,10 

DX< I ) =DX ( I ) +W<  J) #GT  <TRI,J,I) 
CONTINUE 
CONTINUE 

END 


1145  ***t*************************t*****X*************************. 

1146  *  FIND  D<PHI)**2/DX**2 

1147 

1148  SUBROUTINE  PHIXX<TRI ,L1 ,L2,L3,G1 ,G2,G3,DXX) 

1149 

1150  PARAMETER  <MN0DE=151  ,  MNTRIA=50) 

1151  DOUBLE  PRECISION  LI ,L2,L3,G1 ,G2,G3,W< 10) ,DXX< 10) 

1152  DOUBLE  PRECISION  ML<MNTRIA,10,10> 

1153  DOUBLE  PRECISION  NLM<MNTRIA, 16, 10,10) ,NLI (MNTRIA, 4, 10) 

1154  DOUBLE  PRECISION  MG < MNODE , MNODE > 

1155  DOUBLE  PRECISION  GT < MNTRIA, 10, 10) ,LI (MNTRIA, 10 , 4) 

1156  INTEGER  TRI 

1157  COMMON  MG,ML,NLM,NLI,LI,GT 

1158 

1159  U(1)=6.0*L1*G1**2 

1160  14<2)=4'0*L1*G1*G2+2'0*L2*G1**2 

1161  U(3)=4.0*L1*G1#G3+2.0*L3*G1*#2 

1162  U<4)=6.0*L2*G2**2 

1163  14  <  5 ) =4 . 0*L2*G2#G3+2 . 0*L3#G2**2 

1164  U  <  6 ) =4 . 0*L2*G1 *G2+2 . 0*L1 *G2**2 

1165  U(7)=6.0*L3*G3##2 

1166  14  <  8  )  =4 . 0*L3*G  1  *G3+2 . 0*L  1  *G3**2 

1167  14  <  9 )  =4 . 0*L3*G2*G3+2 , 0*L2*G3**2 

1168  W  < 10 ) =2 ♦ 0*  < LI *G2*G3+L2*G1 *G3+L3*G1 *G2 ) 

1169  DO  70  1=1,10 

1170  DXX  < I ) =0  »0 

1171  DO  60  J=1 , 10 

1172  DXX ( I ) =DXX  < I ) +U<  J ) #GT (TRI,J,I) 

1173  60  CONTINUE 

1174  70  CONTINUE 

1175 


1176  END 

1177  *iMtt*********!Mt*#*****##*##*###******#*****#*************** 

1178  *  FIND  D<PHI)**2/CDX*DU) 

1179 

1180  SUBROUTINE  PHIXU<TRI ,L1,L2,L3,F1 ,F2,F3,61 ,G2,G3,DXU) 

1181 
1182 

1183  PARAMETER  <MN0DE=151  ,  MNTRIA-50) 

1184  DOUBLE  PRECISION  F1,F2,F3 

1185  DOUBLE  PRECISION  LI ,L2,L3,G1 ,G2,G3,U< 10) ,DXU< 10) 

1186  DOUBLE  PRECISION  ML(MNTRIA, 10, 10) 

1187  DOUBLE  PRECISION  NLM< MNTRIA, 16, 10, 10) »NLI (MNTRIA, 4, 10) 

1188  DOUBLE  PRECISION  MG ( MNODE  , MNODE ) 

1189  DOUBLE  PRECISION  GT< MNTRIA, 10, 10 > ,LI (MNTRIA, 10,4) 

1190  INTEGER  TRI 

1191  COMMON  MG,ML,NLM,NLI,LI,GT 

1192 

1193  W(1)=6.0*L1*F1*G1 

1194  U(2)=2,0*L1*(F1*G2+F2*G1)+2.0#L2*F1#G1 

1195  U(3)=2,0*L1*(F1*G3+F3#G1)+2»0*L3*F1*G1 

1196  U(4)=6.0*L2*F2*G2 

1197  U  <  5 ) =2 . 0#L2*  (  F3*G2+F2*G3 ) +2 . 0*L3*F2*G2 

1198  U  <  6 ) *2 . 0*L2* ( FI *G2+F2*G1 ) +2 . 0*L1 *F2*G2 

1199  U(7)»6,0*L3*F3*G3 

1 200  U ( 8 ) =2 . 0*L3*  <  F 1 *G3+F3*G1 ) +2 . 0*L1 #F3*G3 

1201  U  <  9 )  =2 . 0*L3*  <  F2*G3+F3*G2 )  +2 . 0*L2*F3*G3 

1202  W< 10 ) =L1#(  F2*G3+F3*G2 ) +L2*<  F1#G3+F3*G1 ) +L3* ( FI *G2+F2*G1 ) 

1203  DO  90  1=1,10 

1204  DXU( I) =0.0 

1205  DO  80  J=l,10 

1206  DXU  < I ) =DXU  < I ) +W  <  J ) #GT (TRI,J,I) 

1207  80  CONTINUE 

1208  90  CONTINUE 

1209 

1210  END 

1211 

1212  *************************************************************** 

1213 

1214  *  DETERMINE  ELEMENT  CASE,  VOLUME,  AND  DERIV'S  OF  TETRAHEDRAL 

1215  *  CO-O  RESPECT  TO  X,  U,  AND  U'  COORDINATES 

1216 

1217  SUBROUTINE  CASEDT( TRI , TRIP , CORDND,PTNODE, TIME ,E ,F ,G,V6, 

1218  C  CASE,U1,U2,U3,X1,X2,X3> 

1219 

1220  PARAMETER  <MN0DE=151  ,  MNTRIA=50 ) 

1221  DOUBLE  PRECISION  ML< MNTRIA, 10 , 10 ) 

1222  DOUBLE  PRECISION  NLM< MNTRIA, 16, 10 , 10 ) ,NLI < MNTRIA, 4 , 10 ) 

1223  DOUBLE  PRECISION  MG < MNODE , MNODE) 

1224  DOUBLE  PRECISION  GT< MNTRIA , 10, 10 ) ,LI < MNTRIA, 10 ,4 ) 

1225  DOUBLE  PRECISION  E< 4) ,F( 4) ,G< 4) ,V6 

1226  DOUBLE  PRECISION  CORDND< MNODE ,2 ) ,D1 , D2 

1227  DOUBLE  PRECISION  U1,U2,U3,B 

1228  DOUBLE  PRECISION  XI ,X2,X3,X2P,U1P,U2P,U3P 

1229  DOUBLE  PRECISION  DE ( 4 , 4) , WK ( 8 ) ,D < 4, 4 ) 

1230  INTEGER  PTNODE< MNTRIA, 1 1 ), TRI , TRIP , CASE , TIME 

1231  COMMON  MG,ML,NLM,NLI,LI,GT 
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1232 

1233 

1234  * 

1235 

1236 

1237 

1238 

1239 

1240 

1241 

1242 

1243 

1244 

1245 

1246  * 

1247 

1248 

1249 

1250 

1251 

1252 

1253 

1254 

1255 

1256 

1257 

1258 

1259 

1260  * 
1261  * 
1262 

1263 

1264 

1265 

1266 

1267 

1268 

1269 

1270 

1271 

1272 

1273 

1274 

1275 

1276 

1277 

1278 

1279 

1280 
1281 
1282 

1283 

1284 

1285 

1286 
1287 


IF  (TIME.EQ.l)  THEN 
CALCULATE  COORDINATES 

X1=C0RDND<PTN0DE<TRI,1> ,1) 
X2=C0RDND(PTN0DE(TRI,4) , 1) 

X3=C0RDND  <  PTNODE <  TRI , 7 ) , 1 ) 
X1P=C0RDND<PTN0DE<TRIP,1>  ,1) 

U1=C0RDND< PTNODE (TRI pi) ,2) 

U2=C0RDND <  PTNODE ( TRI , 4 ) , 2 ) 

U3=C0RDND <  PTNODE  <  TR 1*7) ,2) 
U1P=C0RDND(PTN0DE(TRIP,1) ,2) 

U2P*C0RDND (PTNODE (TRIP, 4) ,2) 

U3P=C0RDND  <  PTNODE ( TR IP , 7 ) , 2 ) 

DETERMINE  THE  CASE  OF  THE  TRIANGLES 
CASE=2 

IF  <X1 .NE.X1P)  THEN 
CASE=1 

IF  (X1.LT.X1P)  THEN 
CAS£=3 
END  IF 

ELSE 

IF  (XI »GT »X2)  THEN 
CASE=>4 
ENDIF 
END  IF 
ENDIF 

ASSEMBLE  THE  COORDINATE  TRANSFORMATION  MATRIX 
ON  CASE 

IF  (CASE.EQ.l)  THEN 
DE(2,1)=X2 
DE(2,2)»X1 
DE(2r3)=*X2 
DE(2,4)=»X1 
DE(3pl)*U2 
DE(3,2)=U1 
DE(3p3)=*U3 
DE(3,4)»U1 
DE(4p 1 )*U1P 
DE(4,2)»U3P 
DE(4p3)*UlP 
DE(4,4)*U2P 
ENDIF 

IF  (CASE.EQ.3)  THEN 
DE(2,1)=X2 
DE<2f2)=X2 
DE(2p3)=Xl 
DE(2,4)=X1 
DE<  3  p  1 )  =>U2 
DE<3p2)=U3 
DE(3p3)=U1 
DE<3p4)=Ul 
DE( 4p 1 ) aUlP 
DE<  4p2)=UlP 
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1288 

DE<4,3)=U3P 

1289 

D£<4,4)=U2P 

1290 

ENDIF 

1291 

1292 

IF  < CASE .EQ .2)  THEN 

“1 

1293 

D£(2,1)=X1 

1294 

DE<2,2)=X2 

; 

1295 

DE(2,3)=X2 

1296 

DE<2,4)=X2 

■. 

1297 

DE(3, 1)=U1 

1298 

DE<3,2)=U3 

■ 

1299 

DE (3,3) =U2 

1300 

DE(3,4)=U3 

1301 

DE(  4, 1 )=U1P 

1302 

DE(4,2)=U3P 

1303 

DE(4,3)=U2P 

j 

1304 

DE(4,4)=U2P 

f  _ 

1305 

IF  (TIME.EQ.2)  THEN 

1306 

DE(3,2)=U2 

1307 

DE(3,3)=U3 

1308 

DE(3,4)=U2 

1309 

DE(4,2)*U2P 

— i 

1310 

DE(4,3)*U3P 

1311 

DE  <  4  *  4 ) *U3P 

— 

1312 

ENDIF 

1313 

ENDIF 

1314  1 

1315 

IF  (CASE.EQ.4)  THEN 

1316 

DE(2,1)=X1 

* 

1317 

DE(2,2)=X2 

V 

1318 

DE(2,3)=X2 

1319 

DE(2,4)=X2 

1320 

DE(3,1)=U1 

, ' 

1321 

DE(3,2)=U3 

1322 

DE(3,3)=U2 

1323 

DE(3»4)=U2 

1324 

DE(4,1)=U1P 

1325 

DE(4,2)=U3P 

1326 

DE(4,3)=U2P 

1327 

DE(4,4)»U3P 

1328 

IF  (TIME.EQ.2)  THEN 

1329 

DE(3,2)=U2 

** 

1330 

DE(3,3)=U3 

1331 

DE(3»4)=U3 

: 

1332 

DE<  4 , 2) =U2P 

1333 

DE(4,3)=U3P 

1334 

DE(4,4)=U2P 

Iz 

1335 

ENDIF 

1336 

ENDIF 

1337 

1338 

DO  10  1=1,4 

- 

1339 

DE< 1 , I) =1 .0 

- 

1340 

10  CONTINUE 

-* 

1341  J 

1342 

*  COPY  MATRIX  TO  AVOID  DECOMPOSITION  BY  IMSL 

i 

1343 

DO  17  1=1,4 

1 


] 
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1344  DO  15  J=1 , 4 

1345  D<I,J>=D£(I,J) 

1346  15  CONTINUE 

1347  17  CONTINUE 

1348 

1349  *  FIND  VOLUME  FROM  MATRIX  DETERMINATE 

1350  Q1=0  *0 

1351  CALL  LINV3F(DE,B,4,4,4,D1,D2,UK,IER) 

1352  V6=D1*2**D2 

1353 

1354  *  DERIVATIVES  OF  NATURAL  COORDINATES 

1355  D1--1 >0 

1356  CALL  LINV3F<D,B,1,4,4,D1,D2,WK,IER) 

1357  DO  20  1=1*4 

1358  E<I)=D<I,2) 

1359  F<I)=D<I,3) 

1360  G<I)=D<I,4) 

1361  20  CONTINUE 

1362 

1363  END 

1364 

1365  ************************************************************************* 

1366  SUBROUTINE  SINFCN<E,F,G, V,SGM,H) 

1367 

1368  PARAMETER  <MN0DE=151  ,  MNTRIA=50> 

1369 

1370  *  FIND  THE  INTERPOLATING  FUNCTION  MATRIX  <20  X  20  FCR  A  CUBIC 

1371  *  IN  3THEN  MULTIPLY  BY  THE  Vi's  TO  GET  THE  Mi's 

1372  *  (THESIS  NOTATION) 

1373  *  RESULT  ARE  THE  BASIS  FUNCTIONS  FOR  THE  TETRAHEDRAL  CUBIC 

1374 

1375 

1376  DOUBLE  PRECISION  ML< MNTRIA, 10, 10 ) 

1377  DOUBLE  PRECISION  NLM (MNTRIA, 16 , 10, 10 ) , NLI < MNTRIA , 4, 10) 

1378  DOUBLE  PRECISION  MG < MNODE , MNODE ) 

1379  DOUBLE  PRECISION  GT< MNTRIA, 10, 10) , LI (MNTRIA, 10,4) 

1380  DOUBLE  PRECISION  M< 17,4,4) ,SGT<20,20) ,UK<8) ,14(4,4) 

1381  DOUBLE  PRECISION  SGM< 5,4 f 4 ) , E < 4 ) , F < 4 ) , G< 4 ) , V< 5 , 20 ) 

1382  DOUBLE  PRECISION  H< 5 , 20 ) , A ,B , C, D1 , D2 ,MT < 4 , 4 ) 

1383  COMMON  MG, ML, NLM, NLI, LI, GT 

1384 

1385  *  ZERO  THE  TETRAHEDRAL  INTERPOLATING  FUNCTION  MATRIX 

1386  DO  10  1=1,20 

1387  DO  5  J=1 , 20 

1388  SGT  < I , J)=0 .0 

1389  5  CONTINUE 

1390  10  CONTINUE 

1391 

1392  *  ASSEMBLE  THE  PARTITIONED  MATRICES  ON  THE  DIAGONAL 

1393  DO  30  K=1 , 4 

1394  DO  20  J»l,4 

1395  M<K,l,J)=l .0*< 1/J) 

1396  M<K,2,J)=3.0*E< J)-<  < J+l >/3)*2.0*E< J) 

1397  M<K,3, J)=3.0*F< J)-< (J+l)/3)*2»  0*F  <  J ) 

1398  M<K,4,J)=3. 0*G <  J ) - < < J+l ) /3 ) *2 . 0*G < J ) 

1399  20  CONTINUE 
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2128  F< 14,5,8)=1 .0 

2129  F<16,5,9>=1.0 

2130  L1=0 .0 


2131 

2132 

2133 

2134 

2135 

2136 

2137 

2138 

2139 

2140 

2141 

2142 

2143 

2144 

2145 

2146 

2147 

2148 

2149 

2150 

2151 

2152  150 

2153 

2154 

2155 

2156 

2157 

2158 

2159 

2160 
2161 
2162 

2163 

2164 

2165 

2166 

2167 

2168 

2169 

2170 

2171 

2172 

2173 

2174 

2175 

2176  160 


L2=l  .0 
L3=0  *0 

CALL  PHIXX<TRI,L1,L2,L3, G1,G2,G3,U1> 

CALL  PHIXU<TRI,L1,L2,L3,F1,F2,F3,G1,G2,G3,U2> 

Ll=1.0 

L2=0  *0 

CALL  PHIXX<TRI,L1,L2,L3,G1,G2,G3,W3) 

CALL  PHIXU<TRI,L1,L2,L3,F1,F2,F3,G1,G2,G3,U4) 
L1=0  .0 
L3=l *0 

CALL  PHIXX<TRI,L1,L2,L3,G1,G2,G3,W5) 

CALL  PHIXU<TRI,L1,L2,L3,F1,F2,F3,G1,G2,G3,W6) 
DO  150  1=1,10 

F(2,I,1) =W3  (I) 

F<3, I , 1 >=W4< I ) 

F<  6,1 , 7>=W5  <  I ) 

F<7,I,7) =W6< I ) 

F  < 10,1 ,4) =U1 < I) 

F(11,I,4) =W2  < I ) 

F< 14,I,7)=W1< I) 

F< 15, I ,7)=W2< I ) 

CONTINUE 
L1=0  »0 
L2=2. 0/3.0 
L3=l. 0/3.0 

CALL  PHIX <  TRI ,L1  ,L2,L3,G1  ,G2 ,G3,I41 ) 

L2=l. 0/3.0 
L3=2. 0/3.0 

CALL  PHI I (TRIP, LI , L2,L3,U2) 

Ll=l. 0/3.0 
L2=2. 0/3.0 
L3-0.0 

CALL  PHIX< TRI ,L1 ,L2 ,L3,G1 ,G2,G3,U3) 

L2=0 .0 
L3=2. 0/3.0 

CALL  PHI I (TRIP,L1 ,L2,L3,U4) 

L3=l. 0/3.0 
L2=l. 0/3.0 

CALL  PHIX < TRI ,L1,L2,L3,G1 ,G2,G3,W5) 

DO  170  1=1.10 

F< 18, I , 10 )=W3< I ) 

F( 20, I , 10) =W5( I ) 

DO  160  J=l,10 

F<17,I,J)=«1(I)*W2( J) 
F<19,I,J)=W5<I)*W4< J) 

CONTINUE 


2177  170  CONTINUE 

2178  ELSE 

2179  IF  (CASE. EQ. 4. AND. TIME. EQ. 2)  THEN 

2180  UU<1)=U1 

2181  UU<  2 ) =U2 

2182  UU ( 3) =U3 

2183  UU  <  4 ) =U3 
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2072 

2073 

2074 

2075 

2076 

2077 

2078 
207? 

2080 

2081 

2082 

2083 

2084  120 

2085 

2086 

2087 

2088 

2089 

2090 

2091 

2092 

2093 

2094 

2095 

2096 

2097 

2098 

2099 

2100 
2101 
2102 

2103 

2104 

2105 

2106  130 

2107 

2108 

2109  140 

2110 
2111 
2112 

2113 

2114 

2115 

2116 

2117 

2118 

2119 

2120 
2121 
2122 

2123 

2124 

2125 

2126 
2127 


L3=l *0 

CALL  PHIXX(TRI,L1,L2,L3,G1,G2,G3,U5) 

CALL  PHIXU<TRI,L1,L2,L3,F1,F2>F3,G1,G2,G3,W6) 
DO  120  1=1,10 
F  <  2 , I , 1 ) =W3  < I ) 

F<3,I,1) =W4< I ) 

F<6, 1, 4-)=Ul<I) 

F  <7,I,4)=W2< I) 

F<10,I,7)=W5<I> 

F<li,I,7)=W6<I> 

F<14,I,7)=U1<I> 

F<15,I,7)=U2<I> 

CONTINUE 
L1=0.0 
L2=2. 0/3.0 
L3=l. 0/3.0 

CALL  PHIX<TRI,L1,L2,L3,G1,G2,G3,W1> 

L2=l, 0/3.0 
L3=2, 0/3.0 

CALL  PHII<TRIP,L1,L2,L3,W2> 

Ll=l. 0/3.0 
L2=0.0 
L3=2. 0/3.0 

CALL  PHII<TRIP,L1,L2,L3,W3> 

L2=2. 0/3.0 
L3=0.0 

CALL  PHIX<TRI,L1,L2,L3,G1 ,  G2,G3,W4) 

L2=l. 0/3.0 
L3=l .0/3.0 

CALL  PHIX <TRI , LI ,L2,L3 ,G1 ,G2,G3,U5 ) 

DO  140  1=1,10 

DO  130  J=l,10 

F<17,I,J)=Ul<r>*W2< J> 

F  < 1 8 , 1 , J ) =W5  < I ) *W3  <  J ) 

CONTINUE 
F(19,I,10) =W4< I ) 

F(20,I,10) =W5  < I ) 

CONTINUE 

ENDIF 

ENDIF 

IF  (CASE.EQ.4.AND.TIME.EQ.1)  THEN 
UU(1)*U1 
UU(2)=U3 
UU  <  3) =U2 
UU<  4) =U2 
F< 1 ,2, 1 >  =  1 ,0 
F<2,2,2)=1.0 
F(4,2,3)=l»0 
F(5,8,7)=1.0 
F<6,8,8)=1 .0 
F<8,8,9)  =  1 .0 
F(?,5,4)=l  .0 
F  < 10,5»5)  =  1 .0 
F< 12,5,6)=!  .0 
F< 13,5, 7) =1 .0 
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2016 

2017 

2018  90 

2019 

2020 
2021 
2022 

2023 

2024 

2025 

2026 

2027 

2028 

2029 

2030 

2031 

2032 

2033 

2034 

2035 

2036 

2037 

2038 

2039 

2040 

2041 

2042  100 

2043  110 

2044 

2045 

2046 

2047 

2048 

2049 

2050 

2051 

2052 

2053 

2054 

2055 

2056 

2057 

2058 

2059 

2060 
2061 
2062 

2063 

2064 

2065 

2066 

2067 

2068 

2069 

2070 

2071 


F(14. I,4)=W5<I> 

F(15,I.4)=W6<I) 

CONTINUE 
L1=0 .0 
L2=l. 0/3.0 
1.3=2 . 0/3.0 

CALL  PHIX(TRI.L1,L2.L3.G1.G2.G3.W1) 

L2=2. 0/3.0 
L3=l. 0/3.0 

CALL  PHI I < TRIP* LI . L2.L3.W2) 

Ll=l. 0/3.0 
L2=2. 0/3.0 
L3=0 . 0 

CALL  PHI  I  <  TRIPOLI  . L2.L3.U3) 

L2=0 .0 
L3=2. 0/3.0 

CALL  PHIX(TRI.L1,L2.L3,G1.G2.G3,W4> 

L3=l. 0/3.0 
L2=l .0/3.0 

CALL  PHIX(TRI .LI . L2.L3.G1 . G2.G3.W5) 

DO  110  1=1^10 

F  < 19 . 1 . 10 ) =W4  < I ) 

F  <  20> 1 . 10 ) =W5  < I ) 

DO  100  J-1.10 

F<17.I.J)=W1 ( I) *U2( J) 

F  < 18. I . J) =W5< I ) *U3( J) 

CONTINUE 

CONTINUE 

ELSE 

IF  ( CASE . EQ  *  2 » AND • TIME .EQ  » 2 )  THEN 
UU(1)=U1 

UU<2)»U2 
UU<3)=U3 
UU<  4)=U2 
F(1.2.1>*1.0 
F<2»2»2)=1 .0 
F(4.2.3)=l *0 

F<  5.5.4 ) =1 . 0 
F<6.5.5)*1 .0 
F(8.5.6)=l .0 
F(9.8.7>*1  .0 
F<10,8.8>=1.0 
F<12.8,9>=1 .0 
F<13,5.7>=1.0 
F< 14. 5. 8) =1.0 
F< 16.5. 9)=1 .0 
L1=0 .0 
L2=l .0 
L3=0 .0 

CALL  PHIXX < TRI »L1.L2»L3.G1»G2»G3»W1) 

CALL  PHIXU<TRI.L1.L2.L3.F1,F2.F3,G1.G2,G3.W2> 

Ll=l .0 

L2=0.0 

CALL  PHIXX(TRI.L1.L2.L3»G1.G2.G3»W3) 

CALL  PH I XU (TRI.L1.L2fL3.fi .F2.F3.G1 .G2.G3.W4) 
L1=0 . 0 
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1980 

1981 

1982 

1983 

1984 

1985 
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1989 
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1992 

1993 
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CALL  PHIX(TRI,L1,L2,L3,G1,G2,G3,U1) 

CALL  PHII<TRIP,L1 ,L2,L3»W6) 

L2=l  *0/3 »0 
L3=0.0 

CALL  PHIX<TRI,L1,L2,L3,G1,G2,G3,W2> 

CALL  PHII<TRIP,L1,L2,L3,U5> 

L1=L2 

L3=L1 

CALL  PHIX<TRI,L1,L2,L3,G1,G2,G3,W3) 

DO  80  1=1,10 

F<  17  *1 ,10 )=W1 < I ) 

F( 18,I,10)=W2( I) 

DO  70  J=l,10 

F<19,I,J)=W3<I)*W5< J> 

F(20,I, J)=W3< I)*W6< J) 

CONTINUE 

CONTINUE 

ENDIF 

IF  < CASE »EQ»2» AND *TIME »EQ » 1 )  THEN 
UU<1)*U1 
UU<2)=U3 
UU(3)*U2 
UU<  4) =U3 
F< 1,2,1)=1 .0 
F<2,2,2>=1.0 
F <  4,2,3>=1 »0 
F<5,8,7)*1.0 
F<6,8,8)=1 »0 
F<8»8*9)=1 *0 
F<  9 ,5 ,4)=1 .0 
F  < 10 ,3,3)  =  1  .0 
F<12,5,6>=1 .0 
F  < 13,8, 4)=1 .0 
F< 14,8,5)=1 .0 
F<16,8,6)=1.0 
L1=0 ,0 
L2=1.0 
L3=0 . 0 

CALL  PHIXX<TRI ,L1 , L2,L3,G1 , G2,G3,W1 ) 

CALL  PHIXU ( TRI ,L1 ,L2,L3,F1 ,F2 ,F3,G1 , G2, G3 , W2) 
Ll=l *0 
L2=0 .0 

CALL  PHI XX (TRI,L1,L2,L3,G1  ,62,63,143) 

CALL  PHIXU <  TRI , LI ,L2, L3,F1 ,F2 ,F3,G1 ,G2,G3 , U4 ) 

L1=0  *0 

L3=1.0 

CALL  PHIXX < TRI ,L1 ,L2 ,L3,G1 ,G2 ,G3,U5) 

CALL  PHIXU ( TRI ,L1 ,L2 ,L3 ,F1 ,F2 ,F3,G1 ,G2 , G3, W6) 
DO  90  1=1,10 

F(2,I,1)=W3(I) 

F<3,I,1) =W4< I ) 

F<  6 , 1 ,7) =W3< I ) 

F<  7 , 1 ,7 ) =U6 < I ) 

F(10,I,4) =W1 ( I ) 

F(11,I,4) =W2  < I ) 

^-39 


h 

1904 

L3=0.0 

1905 

CALL  PHI I < TRIP* LI *L2*L3*U2) 

1906 

L2=0 , 0 

1907 

L3=  1 • 0/3 • 0 

ft 

1908 

CALL  PHII(TRIP*L1*L2,L3*W3> 

IL 

1909 

DO  40  1=1*10 

1910 

DO  30  J=l*10 

■ 

1911 

F < 18* I  * J)=W1 < I )*W2<  J) 

1912 

F(20*I*J)=W1(I)*U3< J) 

1913  30 

CONTINUE 

1914  40 

CONTINUE 

1  , 

1915 

ENDIF 

1916 

“ 

1917 

IF  (CASE *EQ *3)  THEN 

1918 

UU  (  1 )  =U2 

-■* 

1919 

UU ( 2  >  =U3 

^  A 

1920 

UU(3)=U1 

» 

1921 

UU(4)=U1 

.“J 

1922 

F( 1 *5  *  1 )=1 .0 

1923 

F<2*5,2)=1.0 

1924 

F(4*5r3)=l*0 

1925 

F<5»8,1)»1*0 

1926 

F<6*8*2)»1.0 

| 

1927 

F(8,8*3)»1.0 

1928 

F(9*2*7)=l »0 

;■ 

1929 

F  <  10*2,8)=1 *0 

1930 

F< 12*2*9)=1 *0 

1931 

F ( 13*2*4)=1 .0 

M) 

1932 

F  ( 14*2*5) =1 »0 

| 

mJ 

1933 

F(16*2*6)=1.0 

1934 

Ll-0.0 

1935 

L2=l *0 

."-j 

1936 

L3-0.0 

•V 

1937 

CALL  PHIXX<  TRI  *L1  *L2  *  L3  *G1  *  62  *63*141 ) 

1938 

CALL  PHIXU<  TRI *L1 * L2*L3*F1 *F2*F3*G1 *G2  *G3  *  U2) 

i 

1939 

Ll-1.0 

r  - 

•4 

1940 

L2*0  »0 

1941 

CALL  PHIXX  <  TRI *L1 *L2*L3*G1 *G2  *G3* W3) 

1942 

CALL  PHIXU( TRI * LI *L2 *L3*F1 *F2 *F3*G1 *G2 *G3* U4) 

1943 

L1=»0 .0 

1944 

L3-1.0 

1945 

CALL  PHIXX <  TRI *L1 * L2 *L3 , G1 * G2 , G3 , W5 ) 

1946 

CALL  PHIXU(TRI*L1*L2*L3*F1 *F2*F3*G1 , G2*G3*U6) 

■, 

1947 

DO  60  1=1,10 

-  .N 

1948 

F(2*I*1)=U1(I) 

1949 

F<3, I, 1 )=W2< I ) 

-  ■ 

1950 

F(6*I*1)=W5( I) 

1951 

F<  7* I  *  1 )=W6< I ) 

t 

1952 

F(10*I*7)=W3(I) 

-1 

1953 

F< 11  *  I  * 7) =W4< I ) 

1 

1954 

F< 14* I *4) =W3 < I ) 

•j 

1955 

F< 15* I *4) =U4< I ) 

1956  60 

CONTINUE 

'  ■] 

1957 

Ll=2. 0/3.0 

1958 

L2=0 . 0 

1 

1959 

L3=l. 0/3*0 

•  1 
*'  1 

I 

I 


» 
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1848  *  Fi'S  BY  CASE 

184?  IF  ( CASE .EQ »1)  THEN 

1850  UU< 1 ) =U2 

1851  UU<2)=U1 

1852  UU<3)*U3 

1853  UU<4)=U1 

1854  F(l,5,l)*1.0 

1855  F<2,5,2)=1.0 

1856  F<  4,5,3) *1 .0 

1857  F<5,2,7)*1 .0 

1858  F<6, 2,8>*1.0 

1859  F<8,.2,9>*1.0 

1860  F<?,8,1>*1.0 

1861  F< 10»8,2)*1.0 

1862  F(12,8,3>*1.0 

1863  F<13,2,4)*1 .0 

1864  F<14,2,5>=1 .0 

1865  F< 16»2,6)*1 .0 

1866  L1“0 .0 

1867  L2*l  .0 

1868  L3*0 .0 

1869  CALL  PHIXX<TRI,L1,L2,L3,G1,G2,G3,W1> 

1870  CALL  PHIXU<TRI,L1,L2,L3,F1,F2,F3,G1,G2,G3,W2> 

1871  Ll*1.0 

1872  L2-0.0 

1873  CALL  PHIXX<TRI,L1,L2,L3,G1,G2,G3,W3) 

1874  CALL  PHIXU<TRI,L1,L2,L3,F1,F2,F3,G1,G2,G3,W4> 

1875  Ll-0.0 

1876  L3*l *0 

1877  CALL  PHIXX<TRI,L1,L2,L3,G1,G2,G3,U5> 

1878  CALL  PHIXU<TRI,L1,L2,L3,F1,F2,F3,G1,G2,G3,W6> 

1879  DO  10  1*1 t 10 

1880  F<2,I,1 )*W1 <I) 

1881  F<3,I,1)»W2<I> 

1882  F<6,I,7)*W3< I) 

1883  F<7»I»7) *W4  < I > 

1884  F<10,I,1)*W5<I> 

1885  F<11,I,1)=U6<I> 

1886  F( 14, I , 4) *W3<  £ ) 

1887  F ( 15,1 ,4) *W4< I ) 

1888  10  CONTINUE 

188?  Ll*2. 0/3.0 

1890  L2*0»0  . 

1891  L3*l. 0/3.0 

1892  CALL  PHIX<TRI,L1,L2,L3,G1,G2,G3,W1> 

1893  L2»l. 0/3.0 

1894  L3*0 .0 

1895  CALL  PHIX<TRI,L1,L2,L3,G1,G2,G3,W2) 

1896  DO  20  1*1,10 

1897  F< 17, I » 10 )*W1 < I ) 

1898  F(1?,I,10) *W2  < I ) 


1899  20 

CONTINUE 

1900 

Ll*l. 0/3.0 

1901 

L3»l. 0/3.0 

1902 

CALL  PHIX<TRI,L1,L2,L3,G1,G2,G3,W1) 

1903 

Lla2 »0/3*0 
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1816 

1817 

1818 

1819 

1820 
1821 
1822 
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1824 

1825 

1826 

1827 

1828 

1829 

1830 

1831 

1832 

1833 
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*  MULTIPLY  BY  Hi'S  AND  SUM  TO  FIND  THE  SCATTERING  CONTRIBUTION 

*  FROM  THE  <PHI)*<PHI')  TERM 

DO  260  K=l,20 

DO  250  1=1,10 

DO  240  J=1 , 10 

SA<I,J)=SA<I,J)+H<1,K)*F<K,I,J) 

240  CONTINUE 

250  CONTINUE 

260  CONTINUE 


END 

#  SECOND  SCATTERING  INTEGRAL  <  D<PHI)/DX  *  PHI'  ) 

SUBROUTINE  SCATB < U1 ,U2 ,U3, XI , X2,X3,TRI , TRIP , AREAS, H 
C  , CASE, TIME, SB, CORDND,PTNODE) 


PARAMETER  (MN0D£=151  ,  MNTRIA-50) 

DOUBLE  PRECISION  U1 ,U2,U3,X1 ,X2,X3,AREAS<MNTRIA) 

DOUBLE  PRECISION  SB< 10,10) ,U6< 10) 

DOUBLE  PRECISION  A,G1 ,G2,G3,F1 ,F2,F3 
DOUBLE  PRECISION  CORDND  <  MNODE , 2  > 

DOUBLE  PRECISION  Ul< 10) ,W2< 10) ,U3< 10) ,W4< 10) ,W5< 10 ) 
DOUBLE  PRECISION  F<20, 10,10) ,L1,L2,L3,H(5,20) ,UU<4> 
DOUBLE  PRECISION  ML < MNTRI A, 10 , 10 ) 

DOUBLE  PRECISION  NLM<MNTRIA, 16, 10, 10) ,NLI <MNTRIA, 4, 10) 
DOUBLE  PRECISION  MG < MNODE, MNODE) 

DOUBLE  PRECISION  GT(MNTRIA,10,10) ,LI (MNTRIA,10,4) 
INTEGER  CASE, TIME, TRI, TRIP 
INTEGER  PTN0DE<MNTRIA,11) 

COMMON  MG , ML , NLM , NLI , LI ,GT 

*  ZERO  THE  F  MATRICES 

DO  7  K=1 ,20 

DO  6  1=1,10 

DO  5  J»l,10 

F<K,I,J)=0.0 

5  CONTINUE 

6  CONTINUE 

7  CONTINUE 

*  DERIVATIVES  OF  TRIANGULAR  COORDINATES  W.R.T.  SPATIAL 

*  VARIABLES 

A=2 ♦ OfcAREAS  <  TR I ) 

Gl  =  <  U2-U3) /A 
G2*<U3-U1)/A 
G3=<U1-U2)/A 
Fl*<  X3-X2) /A 
F2=<X1-X3)/A 
F3=<X2-X1)/A 

«  ASSIGN  THE  U  COORDS  OF  TETRAHEDRAL  NODES  AND  ASSEMBLE  THE 
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ELSE 


1736 

1737 

1738 

1739 

1740 

1741 

1742  . 

1743 

1744 

1745 

1746 

1747 

1748 

1749 

1750 

1751 

1752 

1753 

1754 

1755 

1756 

1757 

1758 

1759 

1760 

1761 

1762 

1763 

1764 

1765 

1766 

1767 

1768 

1769  180 

1770  190 

1771 

1772 

1773 

1774 

1775 

1776 

1777 

1778 

1779 

1780 

1781  200 

1782 

1783 

1784 


IF  <  CASE . EQ . 4 . AND . TIME . EQ  *  2 )  THEN 

F<1»1>1)-1.0 

F<2fl»2)=1.0 

F<2f2,l>=1.0 

F<3p3p1>*1.0 

F<  4plp3)=l .0 

F(5p4p-4)al  »0 

F<6p5p4>*1 .0 

F<6»4p5)*l *0 

F<7,6r4)*l*0 

F<8,4,6)=1.0 

F<9,7,7>*1 .0 

F<10p8p7)=1.0 

F( 10^ 7f 8)=1 .0 

F<llp9p7)=1.0 

F(12p7p9>»1.0 

F< 13p7p4)  =  l  .0 

F  < 14p8p4) =1 *0 

F<14,7p5>=1.0 

F<15p9p4>-1.0 

F<16p7,6)«1.0 

F<20pl0pl0>=*1.0 

L1=0.0 

L2-1 .0/3.0 

L3=»2. 0/3.0 

CALL  PHII (TRIpLl pL2  pL3p Wl) 
L2«2. 0/3.0 
L3-1. 0/3.0 

CALL  PHII<TRIP,L1pL2pL3pU2> 

DO  190  I»l,10 

DO  180  J-1,10 

F<17pIfJ)=W1<I>*W2( J) 
CONTINUE 
CONTINUE 
Ll*l. 0/3.0 
L2»0.0 
L3=2. 0/3.0 

CALL  PHI I <  TR I  * L 1 » L2  r  L3  p  U2 > 
L2»2» 0/3.0 
L3*0 . 0 

CALL  PHI I (TRIP pLl pL2 pL3pUl ) 

DO  200  I»lplO 

F< 18p I p 10)=W2< I > 
F<19p10pI)»W1<I) 

CONTINUE 

ENDIF 

ENDIF 


1785  *  ZERO  THE  SCATTERING  MATRIX 

1786  DO  230  I»1p10 


1787 

1788 

1789  220 


DO  220  J=1p10 
SA<IpJ)=>0.0 
CONTINUE 


1790  230  CONTINUE 


1680 

CALL  PHII<TRIP,L1,L2,L3,U2> 

1681 

L2=2. 0/3.0 

1682 

L3=0 .0 

1683 

CALL  PHII<TRI,L1,L2,L3,W1> 

1684 

DO  120  1=1,10 

168S 

F<18,10,I)=W2< I) 

1686 

F(19,I,10) =W1( I ) 

1687  120 

CONTINUE 

1688 

END  IF 

1689 

ENDIF 

1690 

1691 

IF  (CASE.EQ.4.AND.TIME.EQ. 1 )THEN 

1692 

F<1,1,1)=1 .0 

1693 

F<2,1,2>=1.0 

1694 

F<2,2,1)=1 .0 

1695 

F<3,3, 1 )=1 *0 

1696 

F<4, 1 ,3>=1 .0 

1697 

F<5,7,7>=1.0 

1698 

F<6,7,8>=1 .0 

1699 

F<6,8,7)=1.0 

1700 

F<7,9,7)=1 »0 

1701 

F<8,7,9>»1,0 

1702 

F<9,4,4)=1 »0 

1703 

F<10,4,5)=1.0 

1704 

F( 10,5,4)*1 .0 

1705 

F<11,6,4)=1.0 

1706 

F  < 12,4,6)=1 »0 

1707 

F(20,10,10)=1.0 

1708 

F<13,4,7)=1 .0 

1709 

F<14,5,7)»1.0 

1710 

F<14,4,8>=1.0 

1711 

F  < 15,6,7)=1 .0 

1712 

F<16,4,9>=1 .0 

1713 

L1=0.0 

1714 

L2=2. 0/3.0 

1715 

L3=l. 0/3.0 

1716 

CALL  PHII(TRI,L1,L2,L3,U1> 

1717 

L2-1. 0/3.0 

1718 

L3=2. 0/3.0 

1719 

CALL  PHII (TRIP,L1 ,L2,L3,W2) 

1720 

DO  160  1=1,10 

1721 

DO  150  J=l,10 

1722 

F<17,I, J)=W1<I)*W2< J> 

1723  150 

CONTINUE 

1724  160 

CONTINUE 

1725 

Ll=l. 0/3.0 

1726 

L2=2. 0/3.0 

1727 

L3=0 .0 

1728 

CALL  PHI I <  TRI ,L1 , L2 , L3, U2) 

1729 

L2=0 .0 

1730 

L3=2. 0/3.0 

1731 

CALL  PHII<TRIP,L1 , L2 ,L3, U1 ) 

1732 

DO  170  1=1,10 

1733 

F< 18,1 , 10 )=U2< I ) 

1734 

F< 19 , 10 , I ) =W1 < I ) 

1735  170 

CONTINUE 

A— 34 

1624 

L3=l. 0/3.0 

1625 

CALL  PHII<TRIP,L1*L2*L3,U2> 

1626 

DO  80  1=1*10 

1627 

DO  70  J=l,10 

1628 

F<17*I*J)=W1<I>*U2< J) 

1629 

70 

CONTINUE 

1630 

80 

CONTINUE 

1631 

Ll=l. 0/3*0 

1632 

L2=2  *0/3*0 

1633 

L3=0»0 

1634 

CALL  PHII  ( TRIP  *L1  *L2  *L3*  142) 

1635 

L2=0 • 0 

1636 

L3=2 *0/3.0 

1637 

CALL  PHII<TRI*L1,L2*L3*W1> 

1638 

DO  90  1=1*10 

1639 

F(18*10»I) =W2< I ) 

1640 

F< 19* I  * 10>=W1 < I ) 

1641 

90 

CONTINUE 

1642 

ELSE 

1643 

IF  < CASE . EQ » 2 « AND • TIME . EQ • 2)  THEN 

1644 

F<1*1*1)»1*0 

1645 

F(2*l ,2)=1 .0 

1646 

F<2*2*1)=1 *0 

1647 

F (3*3*1 )=1 *0 

1648 

F<4, 1 ,3)=1 .0 

1649 

F  <5*4»4)=1 *0 

1650 

F(6*5*4)=l *0 

1651 

F  <6*4»5>=1 *0 

1652 

F<7*6*4)*1 .0 

1653 

F  <8*4*6)=1 *0 

1654 

F<9*7*7)*1 .0 

1655 

F< 10*8*7)=1 .0 

1656 

F<10*7,8)=1.0 

1657 

F<11»9»7)=1*0 

1658 

F< 12,7*9>=1 .0 

1659 

F< 13*4*7)=1 .0 

1660 

F<14*4*8)=1»0 

1661 

F  < 14*5»7)  =  1 ,0 

1662 

F< 1S*6*7)=1 ,0 

1663 

F( 16*4*9)=1 .0 

1664 

F<20*10*10)=1.0 

1665 

L1=0.0 

1666 

L2=2. 0/3.0 

1667 

L3=l. 0/3.0 

1668 

CALL  PHII(TRI*L1*L2*L3*U1) 

1669 

L2=l. 0/3.0 

1670 

L3=2. 0/3.0 

1671 

CALL  PHI I ( TRIP  *L1 , L2  *L3  * W2) 

1672 

DO  110  1=1*10 

1673 

DO  100  J=l*10 

1674 

F<17*I*J>  =W-1  <  I )  *W2  <  J ) 

1675 

100 

CONTINUE 

1676 

110 

CONTINUE 

1677 

Ll-1. 0/3.0 

1678 

L2-0.0 

1679 

L3=2. 0/3.0 

1568 

156? 

1570 

1571 

1572 

1573 

1574 

1575 

1576 

1577 

1578 
157? 

1580 

1581 

1582 

1583 

1584 

1585 

1586  50 

1587 

1588 
158? 
15?0 
15?1 
15?2 
15?3 
15?4  60 
15?5 

_  15?6 

0*  15?7 

15?8 
15?? 
1600 
1601 
1602 

1603 

1604 

1605 

1606 

1607 

1608 
160? 
1610 
1611 
1612 

1613 

1614 

1615 

1616 

1617 

1618 
161? 
1620 
1621 
1622 
1623 


F<?,1,7)*1*0 
F  < 10,1 ,8)=1  .0 
F< 10,2,7)31 *0 
F< 11 ,3,7)=1 .0 
F< 12, 1 ,?) =1 ,0 
F<13,1,4)=1.0 
F<14,1,5>*1.0 
F<14,2,4>*1.0 
F< 15 ,3,4) *1 .0 
F<16,1,6)=1.0 
Ll=2»0/3»0 
L2=0.0 
L3=l. 0/3.0 

CALL  PHI I <TRI .LI , L2,L3,U1 > 
CALL  PHII (TRIPOLI , L2,L3,W2) 

DO  50  1*1,10 

F<17,I,10)=W1<I) 

F<20,10,I)=W2<I) 

CONTINUE 
L2*l. 0/3.0 
L3*0.0 

CALL  PHII<TRI,L1,L2,L3,U1) 
CALL  PHII<TRIP,L1,L2,L3,U2) 

DO  60  1*1,10 

F(18,I,10)*U1<I) 

F<1?,10,I)*W2<I) 

CONTINUE 

ENDIF 

IF  <  CASE • EQ » 2 » AND . TIME • EQ . 1 )  THEN 
F< 1,1,1) *1.0 
F(2,l,2)*1.0 
F  <2,2, 1 )*1 .0 
F (3,3,1 )*1 .0 
F<4, 1 ,3)*1 .0 
F<5,7,7)*1 .0 
F<  6,7 ,8)*1 .0 
F  <6,8,7)*1 .0 
F(7,?,7)*l .0 
F (8,7,?)*1 .0 
F<?,4,4)*1 .0 
F<10,4,5)*1»0 
F(10,5,4)*1.0 
F< 11 , A,4)*l .0 
F<12,4,6)*1.0 
F< 13,7,4)*1 .0 
F< 14,7,5)=1 .0 
F  < 14,8»4)  =  1 .0 
F< 15,?,4)*1 .0 
F<16,7,6)*1.0 
F(20,10,10>*1.0 
Ll-0.0 
L2*l. 0/3.0 
L3»2. 0/3.0 

CALL  PHII<TRI,L1,L2,L3,W1> 
L2*2. 0/3.0 
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1512 

1513 

1514 

1515 

1516 


)  CONTINUE 

)  CONTINUE 

INITIALIZE  THE  Fi'S  DEPENDING  UPON  THE . TETRAHEDRAL  CASE 
AND  TIME 


1517 

IF  (CASE.EQ.l)  THEN 

1518 

F<1,4,1>=1 .0 

1519 

F<2,4,2)=1.0 

1520 

F<2,5,1>=1.0 

1521 

F<3,6,1>*1.0 

1522 

F<4,4,3)=1.0 

1523 

F<5,1,7)=1.0 

1524 

F(6,1,8>“1»0 

1525 

F<6, 2, 71=1.0 

1526 

F<7, 3, 71=1.0 

1527 

F<8, 1,91=1.0 

1528 

F<9,7,11=1 »0 

1529 

F< 10, 8*1) =1*0 

1530 

F<10,7,2)=1 *0 

1531 

F< 11 ,9, 1 )*1 .0 

1532 

F(  12,7,3)=1 *0 

1533 

F<13,1»4)=1.0 

1534 

F<14,1,5)=1*0 

1535 

F<14,2,4)=1.0 

1536 

F( 15,3,4)*1 .0 

1537 

F<16, 1,61*1.0 

1538 

Ll=2. 0/3.0 

1539 

L2-0.0 

1540 

L3=l. 0/3.0 

1541 

CALL  PHIKTRI  »L1  »L2»L3,U1 ) 

1542 

CALL  PHI I (TRIP, LI ,L2,L3,W2) 

1543 

DO  32  1=1,10 

1544 

F(17,I,101=W1(I1 

1545 

F<  20 , 10, 1 )=W2< I ) 

1546  32 

CONTINUE 

1547 

L2=l. 0/3.0 

1548 

L3=0 .0 

1549 

CALL  PHII (TRI , LI , L2,L3,U1 ) 

1550 

CALL  PHII < TRIP, LI ,L2 ,L3,W2> 

1551 

DO  34  1=1,10 

1552 

F(18,10,I)=V»2(I) 

1553 

F< 19,1 , 101=U1 < I ) 

1554  34 

CONTINUE 

1555 

ENDIF 

1556 

1557 

IF  (CASE.EQ.3)  THEN 

1558 

F( 1 ,4, 1 ) =1 .0 

1559 

F<2,4,2)=1 .0 

1560 

F(2,5,l)=l .0 

1561 

F<3,6,1)=1.0 

1562 

F  <  4,4,3)  =  1 .0 

1563 

F<5,7,1 )=1 .0 

1564 

F (6,8, 1 )=1 *0 

1565 

F  <6,7,2)  =  1 .0 

1566 

F(7,9, 1)=1 .0 

1567 

F<8,7,3)  =  1 .0 

1456  SGT<I+12,J+12)=M<12,I,J) 

1457  160  CONTINUE 

1458  170  CONTINUE 


1459 

1460 

1461 

1462 

1463 

1464  180 

1465  190 

1466 

1467 

1468 

1469 

1470  195 


DO  190  1=17,20 

DO  180  J=l,16 

K=<<J-l)/4)+13 
LSJ-(K-13)*4 
SOT (I, J)=M<K,I-16,L) 
CONTINUE 
CONTINUE 

DO  197  1=17,20 

DO  195  J=17,20 

SGT  < I , J)=SGM<5, I-16,J-16) 
CONTINUE 


1471  197  CONTINUE 


1472 

1473  *  FIND  THE  Hi'S 


1474 

1475 

1476 

1477 

1478 

1479 

1480 

1481 

1482 

1483 

1484 

1485 

1486 

1487 

1488 

1489 

1490 

1491 

1492 

1493 

1494 

1495 

1496 

1497 

1498 

1499 

1500 

1501 

1502 

1503 

1504 

1505 

1506 

1507 

1508 

1509 

1510 

1511 


DO  220  1=1,5 

DO  210  J=1 ,20 
H( I , J)=0 *0 
DO  200  L=1 ,20 

H<I,J)»H<I,J)+V<I,L>*SGT<L,J> 

200  CONTINUE 

210  CONTINUE 

220  CONTINUE 

END 

ft*********************************************************************** 

*  CALCULATE  THE  FIRST  <  PHI*PHI'  )  SCATTERING  INTEGRAL 

*  CUBIC  FIT  OVER  A  TETRAHEDRON  WITH  CORNER  FLUXES,  GRADIENTS 

*  AND  FACE  CENTERED  FLUXES  AS  DEGREES  OF  FREEDOM 

SUBROUTINE  SCATA < H ,TRI ,TRIP , CASE , TIME , SA , CORDND , PTNODE ) 

PARAMETER  <MN0DE=151  ,  MNTRIA=50) 

DOUBLE  PRECISION  ML<MNTRIA, 10 , 10) 

DOUBLE  PRECISION  NLM<MNTRIA, 16 , 10 , 10) ,NLI < MNTRIA ,4, 10) 

DOUBLE  PRECISION  MG ( MNODE , MNODE ) 

DOUBLE  PRECISION  GT(MNTRIA,10,10> ,LI(MNTRIA,10,4) 

DOUBLE  PRECISION  X2, X2P,U2,U2P, CORDND < MNODE, 2) 

DOUBLE  PRECISION  F<20, 10, 10) ,U1 < 10) ,W2< 10) ,L1 ,L2,L3 
DOUBLE  PRECISION  SA<10,10) 

DOUBLE  PRECISION  H<5,20) 

INTEGER  TRI, TRIP, CASE, TIME 
INTEGER  PTNODE < MNTRIA, 11) 

COMMON  MG,ML,NLM,NLI,LI,GT 

*  ZERO  THE  Fi ' S 

DO  30  1=1,20 

DO  20  J=l,10 

DO  10  K=1 , 10 

F<I,J,K)=0.0 
10  CONTINUE 
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1400 

A=E<1) 

1401 

B=F< 1) 

1402 

C=G<  1 ) 

1403 

IF  (K.LT.4)  THEN 

1404 

E(1)=E<K+1) 

1405 

F  < 1 ) =F (K+l ) 

1406 

G(1)=G<K+1) 

1407 

E<KF1)=A 

1408 

F<K+1)=B 

1409 

G<KF1)=C 

1410 

ENDIF 

1411 

30 

CONTINUE 

1412 

1413 

*  TAKE 

INVERSES  OF  MATRICES  1-4 

1414 

DO  80  K=l,4 

1415 

DO  50  1=1,4 

1416 

DO  40  J=1 ,4 

1417 

U  <  I ,  J ) =M ( K , I , J ) 

1418 

40 

CONTINUE 

1419 

50 

CONTINUE 

1420 

Dl*-1 »0 

1421 

CALL  LINV3F (U,B»1,4,4,D1,D2,UK,IER) 

1422 

DO  70  1=1,4 

1423 

DO  60  J=1 ,4 

1424 

M<K+8, I ,J)=W< I , J) 

1425 

60 

CONTINUE 

1426 

70 

CONTINUE 

1427 

80 

CONTINUE 

1428 

1429 

*  FIND 

REMAINING  SUB  MATRICES 

1430 

DO  150  K=13, 16 

1431 

DO  110  1=1,4 

1432 

DO  100  J=l,4 

1433 

MT  < I , J ) =0 ♦ 0 

1434 

DO  90  L=1 ,4 

1435 

MT  < I , J)=MT  < I , J)+S6M(K-12, 1 , L) *M(K-4,L , J) 

1436 

90 

CONTINUE 

1437 

100 

CONTINUE 

1438 

110 

CONTINUE 

1439 

DO  140  1=1,4 

1440 

DO  130  J=1 ,4 

1441 

M<K,r,J)=0,0 

1442 

DO  120  L*1 , 4 

1443 

M<K , I , J) =M<  K, I , J)-SGM<  5, 1 ,L) #MT (L , J) 

1444 

120 

CONTINUE 

1445 

130 

CONTINUE 

1446 

140 

CONTINUE 

1447 

150 

CONTINUE 

1448 

1449 

*  ASSEMBLE  INTO  THE  TETRAHEDRAL  (SCATTERING)  INTERPOLATING 

1450 

*  FUNCTION  MATRIX  OF  CONSTANTS 

1451 

DO  170  1=1,4 

1452 

DO  160  J=1 ,4 

1453 

SGT  <I,J)=M<9,I,J) 

1454 

SGT  < 1+4, J+4) =M<10,I,J) 

1455 

SGT ( 1+8, J+8 )=M<11,I,J) 
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2184 

F<1,2,1>»1.0 

2185 

F<2,2,2)»1.0 

2186 

F<4,2,3>=1.0 

2187 

F<5»5»4)=1 .0 

2188 

F<6,5,5>=1.0 

2189 

F<8,5,6)=1.0 

2190 

F<9,8,7>=1.0 

2191 

F< 10»8*8  >»1 .0 

2192 

F<12f8*9>»1.0 

2193 

F< 13*8*4)*1 .0 

2194 

F<14,8,5>=1.0 

2195 

F<16^8,6)=1.0 

2196 

L1*0 .0 

2197 

L2= 1*0 

2198 

L3-0.0 

2199 

CALL  PHIXX<TRI,L1,L2,L3,G1,G2,G3,U1) 

2200 

CALL  PHIXU<TRI,Ll,L2,L3,Fl,F2rF3,Gl,G2,G3,W2> 

2201 

Ll=*l . 0 

2202 

L2=0 .0 

2203 

CALL  PHIXX<TRI,L1,L2,L3,G1,G2,G3,U3> 

2204 

CALL  PHIXU<TRI,Ll,L2,L3rFl,F2,F3,Gl,G2,G3,W4> 

2205 

L1=0.0 

2206 

L3=*l  .0 

2207 

CALL  PHIXX<TRI»L1 fL2.L3fGl .G2fG3. W5) 

2208 

CALL  PHIXU<TRIfLl,L2,L3,Fl,F2,F3rGl,G2>G3,U6> 

2209 

DO  180  Ial t 10 

2210 

F(2pl f l)aW3< I) 

2211 

F<3,lrl)*W4<I) 

2212 

F<6rIr4)»ttl<I) 

2213 

F<7,I,4)»U2<I> 

2214 

F<10,Ir7>»U5<I> 

2215 

F(llrI»7)»W6<I) 

2216 

F<14,I,4>*W5<I> 

2217 

F<15,I,4)a«6<I> 

2218  180 

CONTINUE 

2219 

L1»0.0 

2220 

L2»l. 0/3.0 

2221 

L3»2. 0/3.0 

2222 

CALL  PHIX ( TRIfLl rL2 r L3 pGl . G2. G3.U1 ) 

2223 

L2-2. 0/3.0 

2224 

L3=*l  .0/3 .0 

2225 

CALL  PHI I < TRIP ? LI »L2 . L3  »U2) 

2226 

Ll-l. 0/3.0 

2227 

L2-0.0 

2228 

L3»2. 0/3.0 

2229 

CALL  PHIX<TRI »L1»L2»L3»G1 .G2.G3.U3) 

2230 

L2»2. 0/3.0 

2231 

L3=0.0 

2232 

CALL  PHII<TRIPpLlpL2fL3^W4) 

2233 

L2»l. 0/3.0 

2234 

L3= 1 »  0/3 . 0 

2235 

CALL  PHIX<TRI »L1 »L2»L3 »G1 *G2 » G3» W5) 

2236 

DO  200  I»1 f 10 

2237 

DO  190  J-lrlO 

2238 

F ( 17f I. J) 3W1 ( I ) *U2<  J) 

2239 

F<19,I,J)=W5<I)*U4< J> 
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2240  190  CONTINUE 

2241  F(18,I,10)=U3(I> 

2242  F(20, If 10)=W5(I> 

2243  200  CONTINUE 

2244  ENDIF 

2245  ENDIF 

2246 

2247  *  ZERO  THE  SB  MATRIX 

2248  DO  220  1*1,10 

2249  DO  210  J*l,10 

2250  SB< I , J>=0*0 

2251  210  CONTINUE 

2252  220  CONTINUE 

2253 

2254  *  ASSEMBLE  SB 

2255  DO  250  K=l,20 

2256  DO  240  1*1,10 

2257  DO  230  J=l,10 

2258  SB(I,J)=SB(I,J)+H(2,K)XF(K,I,J)XUU(1)+H(3,K)XF(K,I,J) 

2259  C  XUU(2)+H(4,K)XF(K,I,J)XUU(3)+H(5,K)XF(K,I,J)XUU(4) 

2260  230  CONTINUE 

2261  240  CONTINUE 

2262  250  CONTINUE 

2263 

2264  END 

2265 

2266  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

2267 

2268  SUBROUTINE  PN(PHI ,CORDND, SIGMAS, N,NTRIA, RANGE) 

2269 

2270  PARAMETER  <MN0DE*151  ,  MNTRIA-50) 

2271 

2272  DOUBLE  PRECISION  C0RDND<MN0DE,2) ,PHI (MNODE) ,PSI ( 21 ,9) 

2273  DOUBLE  PRECISION  RANGE, SIGMAS, PERC,TPERC, APERC 

2274  DOUBLE  PRECISION  ML(MNTRIA,10,10) 

2275  DOUBLE  PRECISION  NLM(MNTRIA, 16, 10 , 10) ,NLI (MNTRIA,4, 10) 

2276  DOUBLE  PRECISION  MG ( MNODE, MNODE) 

2277  DOUBLE  PRECISION  GT<MNTRIA, 10, 10 ) ,LI (MNTRIA , 10 ,4) 

2278  INTEGER  N,NTRIA,X,U 

2279 

2280  COMMON  MG,ML,NLM,NLI ,LI ,GT 

2281  X  READ  IN  ARRAY  OF  EXACT  SOLUTION  -  THIS  DATA  FILE  HAS  C*.5 

2282  X  RESULTS 

2283  OPEN  < 18,FILE*'PNDATA5' , STATUS* 'OLD' ) 

2284  REWIND  (18) 

2285  DO  100  1*1,21 

2286  READ  (18,5000)  (PSI( I,U> ,U*1 ,9) 

2287  100  CONTINUE 

2288  CLOSE  (18) 

2289 

2290  X  CALCULATE  PERCENT  DIFFERENCE 

2291  K*0 

2292  TPERC*0  »0 

2293  PRINTX , '  COORDINATES  CURRENTS  FIN  ELEM' 

2294  PRINTX,'  X  U  X  U  FLUX  FLUX  7.  DIFF' 

2295  DO  200  I*1,N-NTRIA,3 


2296 

2297 

2298 

2299 

2300 

2301 

2302 

2303 

2304 

2305 

2306 

2307 

2308 

2309 

2310 

2311 

2312 

2313 

2314 

2315 

2316 

2317 

2318 

2319 

2320 

2321 
EOF, 
EOT. 


IF<INT<C0RDND<I,1)/.25)*.25.EQ.C0RDND<I, 1) >  THEN 
X=NINT <  4*C0RDND  < 1 , 1 ) +1 ) 

U=N INT < 4*C0RDND < I , 2 ) +5 > 

PERC=100.0*ABS<PSI<X,U)-PHI<I>  )/PSI<X,U> 

WRITE  <*,5001)  CORDND< 1,1) ,CORDND( 1,2) ,PHI (I+l),PHI(I+2), 

C  PHI<I)  ,PSI<X,ll)  ,PERC 

IF  < CORDND < 1 , 1). EQ .0 ,0 , AND . CORDND < 1,2) »6E»0 .0) THEN 
GO  TO  200 
ELSE 

IF  <CORDND< 1,1 ) . EQ, RANGE .AND. CORDND < I ,2) »LE. 0.0) THEN 
GO  TO  200 
ENDIF 
END  IF 

TPERC=TPERC+PERC 

K=K+1 

ENDIF 

R ,2399 

200  CONTINUE 

APERC=»TPERC/K 

PRINT*,' AVERAGE  X  DIFFERENCE  IS  ,.',APERC 
D=NTR I A/RANGE 

PRINT*, 'FOR  AN  AVERAGE  OF' ,D, 'TRIANGLES  PER  MEAN  FREE  PATH' 

5000  FORMAT <2X,1P10E13. 4) 

5001  F0RMAT<6< 1X,F7,3) ,2X,F5.2) 

END 


A-45 


181 

===>  CO 

MESHE3.9  MESH 

182 

--»>  XEF9 

183 

IER  IS  . 

.  .  0 

184 

NTRIA 

N  SIGMAS 

185 

46 

151  0.900 

186 

RANGE  IS 

....  3.000000000 

187 

NODAL  VALUES  OF  THE  FLUX 

188 

1 

0.3534  2 

0.3434 

189 

3 

0.3292  4 

0.3044 

190 

5 

0.2500  6 

0.1864 

191 

7 

0.1656  8 

0.1481 

192 

9 

0.1413  10 

0.2554 

193 

11 

0.2288  12 

0.1614 

194 

13 

0.1381  14 

0.1239 

195 

15 

0.0976  16 

0.0904 

196 

17 

0.1745  18 

0.1273 

197 

19 

0.0903  20 

0.0713 

198 

21 

0.0591  22 

0.1208 

199 

23 

0.0999  24 

0.0691 

200 

25 

0.0606  26 

0.0543 

201 

27 

0.0441  28 

0.0398 

202 

29 

0.0826  30 

0 . 0557 

203 

31 

0.0427  32 

0.0377 

204 

33 

0.0337  34 

0.0306 

205 

35 

0 . 0280  36 

0.2539 

206 

ELEMENT 

PENALTY  VALUES 

207 

1 

0.57885E-03 

2 

0.65116E-03 

208 

3 

0.76555E-03 

4 

0.83027E-03 

209 

5 

0.45322E-03 

6 

0.10514E-03 

210 

7 

-.75873E-04 

8 

- .  19484E-03 

211 

9 

-.30669E-03 

10 

-.41093E-03 

212 

11 

- • 731 12E-03 

12 

-.57047E-03 

213 

13 

- .  64482E-03 

14 

- .  48787E-03 

214 

15 

0.39001E-03 

16 

0.51147E-03 

215 

17 

0 . 41353E-03 

18 

0 . 53690E-04 

216 

19 

- , 29261E-04 

20 

- .96223E-04 

217 

21 

- . 21822E-03 

22 

- . 39397E-03 

218 

23 

-.36922E-03 

24 

-.26522E-03 

219 

25 

0.15227E-03 

26 

0.29752E-03 

220 

27 

0. 11491E-03 

28 

0.27448E-04 

221 

29 

-.  16711E-04 

30 

- .33954E-04 

22~> 

31 

- .  12193E-03 

32 

- . 13218E-03 

223 

33 

- . 20394E 

-03  34 

- .87472E-04 

35 

0.90719E 

-04  36 

0 . 10032E-03 

225 

37 

0.64651E 

-04  38 

0 ► 57640E-05 

226 

39 

- . 72660E 

-05  40 

-. 18721E-04 

227 

41 

-. 12605E 

-04  42 

- . 31310E-04 

228 

43 

- .  75572E 

-04  44 

- .31328E-04 

229 

45 

-.  36920E 

-04  46 

- . 48821E-04 

230 

TOTAL  PENALTY  . . . 

.  AND  SUM  OF  ABS< PENALTY) 

ARE  . 

6 

231 

-.46949E 

-04 

0 . 11260E-01 

232 

COORDINATES 

CURRENTS  FIN  ELEM 

233 

X 

U 

X 

U  FLUX 

FLUX 

X  DIFF 

234 

0.000 

1.000 

-0.111 

0.984  0.353 

0.353 

0.00 

235 

0.000 

0.750 

-0.164 

-0.297  0.343 

0.343 

0.00 

236 

0.000 

0.500 

-0.212 

0.115  0.329 

0.329 

0.00 

237 

0.000 

0.250 

-0.373 

0.578  0.304 

0.304 

0.00 

238 

0.000 

0.000 

-0.529 

0.403  0.250 

0.250 

0.00 

239 

0.000 

-0.250 

-0.238 

0.076  0.186 

0.196 

4.67 

240 

0.000 

-0.500 

-0.156 

0.041  0.166 

0.171 

3.06 

241 

0.000 

-0.750 

-0.118 

0.020  0.148 

0.157 

5.42 

242 

0.000 

-1.000 

-0.099 

-0.161  0.141 

0.147 

3.56 

243 

0.750 

1.000 

-0.119 

0.552  0.255 

0.262 

2.61 

244 

0.750 

0.750 

-0.141 

0.045  0.229 

0.236 

3.12 

245 

0.750 

0.250 

-0.129 

0.148  0.161 

0.169 

4.49 

246 

0.750 

0.000 

-0.234 

0.193  0.138 

0.142 

2.79 

247 

0.750 

-0.250 

-0.129 

0.133  0.124 

0.125 

0.99 

248 

0.750 

-0.750 

-0.073 

0.040  0.098 

0.101 

3.28 

249 

0.750 

-1.000 

-0.057 

-0.064  0.090 

0.091 

1.23 

250 

1.500 

1.000 

-0.079 

0.173  0.174 

0.184 

5.08 

251 

1.500 

0.500 

-0.082 

0.158  0.127 

0.132 

3.38 

252 

1.500 

0.000 

-0.116 

0.090  0.090 

0.095 

4.52 

253 

1.500 

-0.500 

-0.058 

0.043  0.071 

0.075 

4.52 

254 

1.500 

-1.000 

-0.038 

0.020  0.059 

0.062 

4.35 

255 

2.250 

1.000 

-0 . 060 

0.122  0.121 

0.127 

4.74 

256 

2.250 

0.750 

-0.054 

0.063  0.100 

0.106 

5.62 

257 

2.250 

0.250 

-0.047 

0.044  0,069 

0.073 

5.56 

258 

2.250 

0.000 

-0.094 

0.075  0.061 

0.063 

4.46 

259 

2.250 

-0.250 

-0.052 

0 . 066  0 . 054 

0.056 

3.07 

260 

2.250 

-0.750 

-0.030 

0.032  0.044 

0.045 

3.03 

261 

2.250 

-1.000 

-0.022 

0.039  0.040 

0.042 

4.30 

262 

3.000 

1.000 

-0.040 

0.076  0.083 

0.087 

4.85 

263 

3.000 

0.500 

-0.030 

0.045  0.056 

0.058 

4.28 

264 

3.000 

0.000 

— 0 . 035 

0.012  0.043 

0.043 

0.00 

265 

3.000 

-0.250 

-0.033 

0.291  0.038 

0.038 

0.00 

266 

3.000 

-0.500 

-0.018 

-0.015  0.034 

0.034 

0.00 

267 

3.000 

-0.750 

-0.017 

0.045  0.031 

0.031 

0.00 

268 

3.000 

-1.000 

-0.013 

0.070  0.028 

0.028 

0.00 

269 

AVERAGE 

X  DIFFERENCE  IS  . 

.  -  3.877686025 

270  FOR  AN  AVERAGE  OF  15.333  TRIANGLES  PER  MEAN  FREE  PATH 


ED 

C90UT 

LI, 

1 ,300 

1 

===>  CO 

MSHE3.9C  MESH 

r> 

===>  XE 

3 

IER  IS  . 

.  .  0 

4 

NTRIA 

N  SIGMAS 

5 

46 

151  0.900 

6 

RANGE  IS 

»  »  »  »  3  #000000000 

7 

NODAL  VALUES  OF  THE  FLUX 

8 

1 

0 . 3534  2 

0.3434 

9 

3 

0.3292  4 

0.3044 

10 

5 

0.2500  6 

0.1892 

11 

7 

0.1691  8 

0.1519 

12 

9 

0.1439  10 

0.2580 

13 

11 

0.2325  12 

0.1660 

14 

13 

0.1419  14 

0.1269 

15 

15 

0.1004  16 

0.0927 

16 

17 

0.1818  18 

0.1293 

17 

19 

0.0930  20 

0.0736 

18 

21 

0.0604  22 

0.1253 

19 

23 

0.1035  24 

0.0714 

20 

25 

0.0626  26 

0 . 0562 

21 

27 

0.0450  28 

0.0406 

22 

29 

0.0856  30 

0.0571 

23 

31 

0.0427  32 

0.0377 

24 

33 

0.0337  34 

0.0306 

25 

35 

0.0280  36 

0  •  2804* 

26 

ELEMENT 

PENALTY  VALUES 

27 

1 

0.69119E-03 

2 

0.79918E-03 

28 

3 

0.70362E-03 

4 

0.79296E-03 

29 

5 

0.40508E-03 

6 

0.10532E-03 

30 

7 

-.84340E-04 

8 

- ♦ 19188E-03 

31 

9 

-.30460E-03 

10 

-.42177E-03 

32 

11 

- .75959E-03 

12 

-.58829E-03 

33 

13 

-.66308E-03 

14 

- . 51004E-03 

34 

15 

0.44435E-03 

16 

0 .54365E-03 

35 

17 

0 . 40379E-03 

18 

0 . 57488E-04 

36 

19 

-.32378E-04 

20 

- • 10216E-03 

37 

21 

- . 23084E-03 

22 

- . 41572E-03 

38 

23 

- . 39299E-03 

24 

-.28130E-03 

39 

25 

0 . 16749E-03 

26 

0 . 31947E-03 

40 

27 

0  » 1 1766E-03 

28 

0 . 30371E-04 

41 

29 

- . 17476E-04 

30 

- . 3601 IE— 04 

42 

31 

- . 12910E-03 

32 

- . 14093E-03 

43 

33 

- .  22000E-03 

34 

-.951 10E-04 

44 

35 

0.99120E-04 

36 

0 . 10831E-03 

45 

37 

0.67810E-04 

38 

0 , 71457E-05 

46 

39 

- . 77406E-05 

40 

- • 20200E-04 

47 

41 

-. 18113E-04 

42 

- . 26107E-04 

48 

43 

-.80449E-04 

44 

- . 35963E-04 

49 

45 

-.40512E-04 

46 

- . 54356E-04 

50 

TOTAL  PENALTY  ....  AND 

SUM 

1  OF  ABS< PENALTY)  ARE  ►. 

51 

-.39050E-04 

0 . 11767E-01 

52 

COORDINATES  CURRENTS  FIN  ELEM 

53 

X 

U  X 

U  FLUX  FLUX  Z  DIFF 

54 

0.000 

1.000  -0.134 

0.040  0.353  0.353  0.00 

Ap  p  end i x  B  . 

The  Absorbing  Term  (quadratic,  and  cubic) 
The  absorbing  term  is 


(3-2) 


Uhich  may  be  written  as 


<£  M  A  6 1 


(3-5) 


I  e 


B 


The  matrices  on  subsequent  pages  Labeled  as  QCOABS  and 

CCOABS  (for  Quadratic  Coefficients  Absorbing  term  and  Cubic 

Coefficients)  generate  MA  of  equation  (3-5)  from 

■u 


M*  -  QcofiiSS 


(  B  - 1 ) 


and  in  the  cubic  case  from 


/MA  ~ 


-  Z4 


CCOfl-TSS 


?! 


(  B  -  2  ) 


where  A  is  triangle  area. 


3-1 


.  .-i.~Z.M-m  -1-«  Xm  -1  « 


cat  coabs 


QCOABS 


24.0 

4.0 

4.0 

6.0 

2.0 

6.0 

4.0 

24.0 

4.0 

6.0 

6.0 

2.0 

4.0 

4.0 

24.0 

2.0 

6.0 

6.0 

6.0 

6.0 

2.0 

4.0 

2.0 

2.0 

2.0 

6.0 

6.0 

2.0 

4.0 

2.0 

6.0 

2.0 

6.0 

2.0 

2.0 

4.0 

cat  ccaabs 


D« 


CCOABS 


720.0 

120.0 

120.0 

36.0 

12.0 

120.0 

48.0 

24.0 

48.0 

12.0 

120.0 

24.0 

48.0 

12.0 

8.0 

36.0 

48.0 

12.0 

720.0 

120.0 

12.0 

12.0 

8.0 

120.0 

48.0 

48.0 

36.0 

12.0 

120.0 

24.0 

36.0 

12.0 

48.0 

36.0 

48.0 

48.0 

12.0 

36.0 

12.0 

12.0 

12.0 

8.0 

12.0 

48.0 

36.0 

24.0 

12.0 

12.0 

24.0 

12.0 

48.0 

36.0 

48.0 

12.0 

24.0 

36.0 

12.0 

12.0 

8.0 

12.0 

12.0 

48.0 

36.0 

12.0 

12.0 

120.0 

36.0 

12.0 

48.0 

24.0 

24.0 

48.0 

12.0 

36.0 

12.0 

48.0 

12.0 

8.0 

12.0 

12.0 

12.0 

720.0 

120.0 

120.0 

24.0 

8.0 

120.0 

48.0 

24.0 

12.0 

12.0 

120.0 

24.0 

48.0 

12.0 

12.0 

24.0 

12.0 

12.0 

8.0 

3-2 


ID< 


-  -V  -y  ’ 


Appendix  C  -  Boundary  Term  (quadratic  and  cubic) 

The  boundary  term  (3-6)  is 

u.  ^4  r  t_  (^4  ^  QT  6T^ 

J  '  ~  “  (3-6) 


When  u  is  expanded  as  in  equation  2-3,  the  boundary  matrix 
can  be  written  as  the  sum  of  3  matrices. 


-  _L  (£ 


A(Bl  z'42.2-  +•  3 


6T  ^ 


(3-10) 


^  r 

— y 


1 


These  matrices  are  complicated  by  the  fact  that  the  basis 
functions  contain  derivatives  of  natural  co-ordinates,  which 
are  distinct  for  every  separate  triangle  geometry.  In  the 
quadratic  caae/wy  Is  given  bv 

2^3*- 

A  3- 

L  A  3  3  *•  3 '  J 

-A— 

the  product  /ni  for  any  of  the  three  matrices  results  in  what 
can  be  thought  of  as  a  matrix  of  constants,  multiplied  by  3  J  S 
as  appropriate.  If  a  matrix  ^  i3  formed  of  the  <3  c  $  . 


(2-19) 


b  - 


3, 

0 

a 

33 

O 

3' 

3  2. 

-'X 

33 

a  T 

3' 

(  C  -  1 ) 


'J 


*  i 


.i 

h 


U 


> 


x 


overlayed 


then  this  can  be  computed  for  each  triangle  and  "  overlayed  " 
in  a  seise  on  each  column  of  constants  to  produce  the  boundary 
matrix.  As  an  example  consider  the  matrix  referred  to  on  the 
next  page  as  QCOBNDI  (  for  Quadratic  COeffients  Bndry  term  #1  ) 
It  is  a  6  X  12  matrix,  that  when  D  is  overlayed  on,  and  when 

multiplied  by  J4,  2.  A  produces  MBI 

(a  l  f  ,  fO 

^3i  *  ^33 

* 

MB 2  and  MB3  must  be  farmed  uf  course  with  the  appropriate 
constants  ( ^  and  ^  )  from  the  matrices  labeled  QC0BND2  and 
QC0BND3 . 


MTU  - 


M,  f*2  A 


<o 


(  C  -  2  ) 


Boundary  matrices  for  the  cubic  fit  are  found  in  an 
analagous  manner,  with  3  (10  x  20  )  matrices,  CC0BND1  ,  CC0BND2, 
and  CC0BND3.  In  the  cubic  case,  the  derivative  matrix  D  to  be 


over layed  i s 


O 

9/ 

3. 

93 

(3 

33 

3' 

3j 

3 

3j 

3. 

u3 

32 

J1 

(  C  -  3  ) 


la  this  case  row  10  must  be  augmented  by  another  term.  T 


Last  row  of  CCOBNDl,  CCOBND2 ,  and  CC0BND3  represent  3  each 
dimension  (10)  matrices  (BR1,  BR2,  BR3).  After  is  formed 

from 


MJS  —  MJZ\  1-  2.  +■  A1  BS 


(C-4 


if 

For  i=l , 10 

=  r$tf  \  C.i)f>r  W,  v  +■  f  HftZCc)  ^  ^  3 

Then 

For  1-1,10 

nrs0otl)  =  A *2(10,  i)  +  HR  c0  *  ‘2-'4  *  6  ' 


and  the  boundary  matrix  is  now  completely  formed. 


cat  cobndl 


48.0 

.0 

8.0 

12.0 

.0 

12.0 

12.0 

.0  ' 

4.0 

6.0 

24.0 

6.0 

6.0 

6.0 

2.0 

24.0 

6.0 

4.0 

7.  cat 

cobnd2 

12.0 

.0 

12.0 

8.0 

.0 

48.0 

4.0 

.0 

12.0 

4.0 

6.0 

24.0 

2.0 

4.0 

6.0 

6.0 

2.0 

6.0 

7.  cat  cobnd3 

QC0BND1 

,0  8.0  .0 

.0  4.0  .0 

,0  12.0  .0 

4.0  2.0  4.0 

6.0  6.0  2.0 

2.0  4.0  6.0 


QC0BND2 

.0 

4.0 

.0 

.0 

8.0 

.0 

.0 

12.0 

.0 

6.0 

4.0 

2.0 

24.0 

6.0 

4.0 

6.0 

2.0 

6.0 

12.0 

.0 

4.0 

8.0 

.0 

4.0 

4.0 

.0 

4.0 

4.0 

6.0 

2.0 

2.0 

4.0 

2.0 

6.0 

2.0 

2.0 

8.0 

.0 

4.0 

12.0 

.0 

12.0 

4.0 

.0 

8.0 

6.0 

4.0 

6.0 

2.0 

6.0 

4.0 

4.0 

2.0 

2.0 

. 0  12.0  .0 

.0  4.0  .0 

.0  8.0  .0 

2.0  2.0  6.0 

2.0  4.0  2.0 

2.0  6.0  4.0 


.0  4.0  .0 

.0  4.0  .0 

.0  4.0  .0 

2.0  2.0  2.0 

6.0  2.0  2.0 

4.0  2.0  2.0 


0C0BND3 


12.0  .0  4.0 

4.0  .0  12.0 

8.0  .0  8.0 

2.0  6.0  6.0 

4.0  2.0  4,0 

6.0  4.0  2.0 


.0 

12.0 

.0 

.0 

12.0 

.0 

.0 

48.0 

.0 

2.0 

6.0 

6.0 

6*0 

24.0 

6.0 

4  ♦  0 

6.0 

24.0 

4.0  .0  4.0 
4.0  .0  8.0 
4.0  .0  12.0 
2.0  2.0  4.0 
2.0  2.0  6.0 
2.0  2.0  2.0 


.0  8.0  .0 

,0  4.0  .0 

.0  12.0  .0 

2.0  2.0  4.0 

4.0  6.0  2.0 

6.0  4.0  6.0 


G-4 


7.  cat  ccostr3 


CC0STR3 


432.0 

.0 

.0 

72.0 

.0 

144.0 

.0 

216.0 

.0 

144.0 

.0 

72.0 

.0 

.0 

72.0 

.0 

24.0 

.0 

72.0 

.0 

24.0 

.0 

432.0 

.0 

.0 

216.0 

.0 

144.0 

.0 

72.0 

.0 

144.0 

.0 

72.0 

144.0 

.0 

32.0 

24.0 

24.0 

48.0 

48.0 

72.0 

74 . ") 

48.0 

72.0 

24.0 

.0 

48.0 

24.0 

24.0 

8.0 

32.0 

...  -*  ♦  0 

24.0 

8.0 

144.0 

144.0 

.0 

48.0 

72.0 

48.0 

48.0 

48.0 

24.0 

48.0 

48.0 

216.0 

144.0 

.0 

48.0 

24.0 

72.0 

48.0 

192.0 

72.0 

72.0 

48.0 

72.0 

24.0 

.0 

96.0 

24.0 

24.0 

8.0 

48.0 

24.0 

24.0 

8.0 

720.0 

144.0 

.0 

192.0 

72.0 

240.0 

48.0 

96.0 

24.0 

240.0 

48.0 

72.0 

.0 

.0 

72.0 

.0 

24.0 

.0 

72.0 

.0 

24.0 

.0 

432.0 

.0 

.0 

216.0 
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with  groups  of  2  or  3  columns  at  a  time,  as  divided  by  dotted 
Lines  on  fig  C-l,  until  all  10  columns  of  the  streaming  matrix 
are  assembled. 

In  the  cubic  case,  the  arrays  listed  as  CCOSTRI,  CC0STR2, 

CC0STR6  represent  10  x  33  matrices  of  constants,  for  rows  1 
through  9  of  the  6  streaming  matrices  with  a  (1  X  18)  row 
matrix  below  to  augment  row  10  terms.  DS  is  a  10  x  33  matrix 
of  ,  which  when  overlayed  on  the  sum  of  CCOSTRI  thru  6, 

multiplied  by  the  appropriate  U  'S  and  (  )  from  the  integra¬ 

tion,  forms  the  streaming  matrix. 

Column  1-3  of  DS  overlayed  on  columns  1-3  of  the  COSTR  sum 
produce  the  first  column  of  MS.  Subsequent  columns  of  the 
streaming  matrix  are  founed  by  the  next  3  or  4  columns  of  DS, 
overlayed  on  the  corresponding  CCOSTR  columns,  as  separated  by 
the  dotted  lines  in  DS  of  figure  c-2. 

Row  10  of  the  streaming  matrices  must  be  augmented  by  the 

dimension  (18)  matrix  ( SRI , SR2 , . .  ,  SR6 )  on  the  Last  two  lines  of 

COSTR,  such  that 
F  or  t:i  f  t  % 

s*cc)z f^u^s^cc)  s/nc-) tic :}  +  3^^  S#9CC) 

"  2A/9z 
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msc,o,i)+  + 

z  *sC/o, 3J  -fcsve^J  5/ 

MS  C'0,Sj  Z  /tSCmS-) 
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S)  -  A*s(jo,%)  -rS^Jq,2-  S*0o)  1,  3  ! 

MSt’O'Q)  r  M  SO  0,1)  *•  S%02) 

/*  sC'°,'a)  -  ~  i,1 4.  5  zL  '-'jcj'O  z  *  C:  3,  3  3 


2. 


Appendix  D  -  The  Streaming  Term  (quadratic  and  cubic) 


The  Streaming  term  is 


(D-l) 


which  can  be  written  as  the  sum  of  6  distinct  matrices 

—  -L  (pJ-  f/uS  i  A^S-T +MS&T]  ^ 

-  —  £  £r  /i^  ^_r  (3-13) 

Evaluating  these  mtrices  involves  taking  the  product 

which  results  in  cross  products  of  the  natural  coordinates 

derivatives  with  respect  to  cartesian  coordinates.  The  six 

streaming  matrices,  as  in  the  boundary  case  (Appendix  C)  can  be 

thought  of  as  distinct  matrices  of  constants,  which  after  being 

7.  1 

multiplied  respectively  by  M,  /  f  /U.-$  / 

and  2/J,jaz  can  be  summed,  and  then  "overlayed"  by  a  matrix  of 
.  Due  to  the  cross  products,  this  matrix  of  derivatives 
( PS )  is  complicated.  It  Is  generated  in  Subroutine  Stream, 
(appendix  A)  for  the  cubic  case,  and  written  out  below  for  both 
the  quadratic  and  cubic  cases. 

DS  for  the  quadratic  case  is  listed  in  figure  d-l.  After 
multiplying  the  6  QCOSTR  matrices  by  the  appropriate  u  values; 
and  a  factor  of  <  )  from  the  integration,  they  may  be 

summed  to  form  a  single  6x18  matrix.  The  first  two  columns  of 
this  matrix  overlay  on  the  first  two  columns  of  DS  to  form  the 

first  column  of  the  streaming  matrix.  The  process  continues 

D-1 
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Appendix  E  -  The  Scattering  Integrals 
The  first  scattering  integral  is 


</u  J  u/  DC  (p  & 
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where  PC  *  E-Z-  ~  ^5  .  If  the  product  (ft  —  F 

then  a  cubic  can  be  specified  over  the  tetrahedron  using  the 


twenty  degrees  of  freedom  specified  in  figure  2-6.  Consider 
case  1  as  depicted  in  figure  4-4. 
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3  '  1  3  3  ’ 

these  are  not  finite  element  interpolation  nodes,  but  that  they 
can  be  written  in  terms  of  these  nodes  using  (2-13).  The  second 
scattering  integral  is 
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Continuing  to  number  as  in  figure  4-4,  the  integrals  for 
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For  case  4,  points  11,12,15  and 

(-L  ol)  (0,1  x  )  and  (0,0-  2). 

s  3  3  3  '  3 1  3 

(^,£  O  )  (-l-oJ;),  (oi.J_)  and 

ij1  3'  '3  ' 3,3 

triangle . 


16  are  local  at  (?,i  0  ) 

3  '  3  ' 

Points  13,14,17  and  18  are  at 
(  O, -W  —  )  on  the  non  local 

i  3>  $ 
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Finite  Element  Meshes 


1.  Streaming  meshes 

Meshl 
Mesh2 
Mesh3 
Mesh4 
Mesh5 
Mes  h6 

2.  Scattering  meshes 

Mesha 

Meshb 

Meshc 

Meshd 

Meshes  1  through  4  are  identical  to  those  of  reference  (2). 
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34 

4 

53 

54 

46 

35 

27 

28 

4 

55 

COLUMN 

FIRST  ELEMENT  NUMBER 

OF 

56 

1 

1 

14 

57 

2 

15 

10 

58 

3 

25 

10 

59 

60 

4 

35 

12 

61 

NODE 

X-AXIS  U-AXIS 

62 

1 

0.000  1.000 

63 

2 

0.000  0.750 

64 

3 

0.000  0.500 

65 

4 

0.000  0.250 

66 

5 

0.000  0.000 

67 

6 

0.000  -0.250 

68 

7 

0.000  -0.500 

69 

8 

0.000  -0.750 

70 

9 

0.000  -1.000 

71 

10 

0.250  1.000 

72 

11 

0.250  0.750 

73 

12 

0.250  0.250 

74 

13 

0.250  0.000 

75 

14 

0.250  -0.250 

76 

15 

0.250  -0.750 

77 

16 

0.250  -1.000 

78 

17 

0.500  1.000 

79 

18 

0.500  0.500 

80 

19 

0.500  0.000 

81 

20 

0.500  -0.500 

82 

21 

0.500  -1.000 

83 

22 

0.750  1.000 

84 

23 

0.750  0.750 

85 

24 

0.750  0.250 

86 

25 

0.750  0.000 

87 

26 

0.750  -0.250 

88 

27 

0.750  -0.750 

89 

28 

0.750  -1.000 

90 

29 

1 . 000  1 . 000 

91 

30 

1.000  0.500 

92 

31 

1.000  0.000 

93 

32 

1.000  -0.250 

94 

33 

1.000  -0.500 

95 

34 

1.000  -0.750 

96 

97 

35 

1.000  -1.000 

98 

99 

NODE  < NB) 

20 

FLUX 

100 

1 

1.0142E+00 

101 

3 

9.8832E-01 

102 

4 

7 . 6712E-01 

103 

6 

9.7586E-01 

104 

7 

5.2627E-01 

105 

9 

9.8330E-01 

106 

10 

2.7547E-01 

107 

12 

8 . 1002E-01 

108 

13 

1 . 2126E-01 

109 

15 

6 . 1684E-01 

110 

91 

6 . 5917E-03 

?-26 

c« 


Appendix  G  -  Subroutines  to  numerically  evaluate  the  scattering 
integral 

Glossary  of  Variables 

xl , x2 , x3 , ul ,  u2  ,  u3  -  triangle  geometric  coordinates 
x( 49 )  ,  u ( 4 9  )  -  coordinates  of  49  integration  points 

MC  -  coefficient  matrix  of  eqn  2-4 
DET  -  deteminate  of  MC 
DU(7)  -  delta  u  at  each  of  the  seven 
spatial  integration  points 
LX(3)  -  derivatives  of  -&i f  and  UJ ,  Y,  “C  >  ^ 

DX(mntria)  -  column  width 

L  (49,3)  -  array  storing  natural  coordinates  of  integration 
points 

M ( 4 9 , 1 0 )  -  array  storing  m  of  2-29  evaluated  at  the  integration 


FLUX(49,10)  -  flux  at  integration  points 
DFLUX( 49 , 10 )  -  derivative  of  flux 


U 1 1 , U 1 2 , .  . .  , U 1 7  -  u  integrals  at  the  seven  points  needed  for 


spatial  integration 


UtJHO'«fflMO'«>WHHO'Oa)\IO'«*UMH‘0'Ofl)''JO‘UUWMHO'OOOMO"OIJ>WWHO 


LI, 1 ,160 

1  **########*##:M##########**#*##*##**##******#####****##******* 

3 

4  SUBROUTINE  LCORD ( AREAS »TRI , PTNODE ,CORDND ,M ,MX , DU,DX,U,X) 

5 

6  PARAMETER  <MN0DE=151  ,  MNTRIA=46) 

7  DOUBLE  PRECISION  XI ,X2, X3 , U1 , U2 ,U3 , XX ( 7) , X ( 49 ) ,U< 4?) 

8  DOUBLE  PRECISION  MC<3,3) ,DET,DU<7) ,LX<3) ,DX(MNTRIA) ,L(49,3) 

9  DOUBLE  PRECISION  AREAS< MNTRIA) ,CORDND<MNODE, 2) , M( 49 , 10 ) ,MX ( 49 , 
DOUBLE  PRECISION  ML ( MNTRIA, 10, 10 ) ,MG ( MNODE ,MNODE) 

DOUBLE  PRECISION  NLM< MNTRIA, 14, 10 , 10 ) ,LI < MNTRI A, 10 , 7) 

DOUBLE  PRECISION  NLI <MNTRIA, 7 , 10 ) , AS ( MNODE* < MNODE- 1 ) /2) ,F ,G 
INTEGER  CASE,TRI ,PTNODE( MNTRIA, 11 ) 

COMMON  MG, ML, NLM, LI, NLI, AS 

*  THIS  SUBROUTINE  FINDS  THE  NATURAL  COORDINATES  NEEDED  FOR 

*  NUMERICAL  INTEGRATION  OF  THE  SCATTERING  INTEGRAL 

#  49  POINTS  FOR  WEDDLES  N-6 

#  GET  THE  <X,U)  COORDINATES  OF  THE  TRIANGLE 
X1=C0RDND< PTNODE  ( TRI  ,  1  >  ,  1  > 

X2-C0RDND  <  PTNODE ( TRI , 4  > , 1 > 

X3-C0RDND  <  PTNODE ( TR 1 , 7  )  ,  1 ) 

U1=C0RDND<PTN0DE(TRI , 1 > ,2) 

U2-C0RDND ( PTNODE  <  TRI , 4 ) , 2 ) 

U3=C0RDND< PTNODE (TRI ,7) ,2) 

*  DETERMINE  THE  ORIENTATION  OF  THE  ELEMENT 
CASE-2 

IF  <X1,GT.X2>  THEN 
CASE-1 
END  IF 

#  GET  X  COORDS  OF  NUMERICAL  INTEGRATION  POINTS 
XX < 1 )-MIN<Xl , X2,X3) 

XX<7)-MAX<X1 , X2,X3) 

F— <XX<7)-XX< 1 ) )/6*0  * 

DO  50  1=2,6 

XX(I)=XX< 1)+(I-1 >*F 
50  CONTINUE 

DO  70  1=1,7 
J=7*I-6 
DO  60  K-J , J+6 
X  <  K ) =XX  < I ) 

60  CONTINUE 

70  CONTINUE 


*  GET  U  COORDS  OF  THE  SAME  POINTS 
IF  <  CASE • EG • 1 )  THEN 
U< 1 )-U3 
U  <  7 ) =U2 
F=  <  U1-U3 ) /6  *0 

54  G=  <  U1-U2 ) /6 .0 

55  DO  80  1=8,36,7 


G-2 


36 

J=(I-l>/7 

a 

57 

U<I)=U3+J#F 

58 

U< I+6)=U2+J*G 

59  80 

CONTINUE 

60 

U ( 43 ) =U1 

61 

U<49)=U1 

' 

62 

ELSE 

63 

U< 1)=U1 

• 

64 

U<7)=U1 

_ 

65 

F=<U2-Ul)/6.0 

66 

G=<U3-Ul)/6*0 

67 

DO  95  1=8,36,7 

68 

J=<I-l>/7 

69 

U<I)=U1+J*F 

70 

U< I+6)=U1+J*G 

71  95 

CONTINUE 

( 

72 

U  <>  43  )  =U2 

-1 

73 

U<49)=U3 

74 

END  IF 

75 

DO  100  1*1,43,7 

76 

F=(U< I+6)-U< I) )/6.0 

77 

DO  98  J*1 ,5 

78 

U<I+J>=U<I)-KJ*F 

79  98 

CONTINUE 

80  100 

CONTINUE 

31 

*, 

82  *  COMPUTE  THE  LOCAL  NATURAL  COORDINATES 

jj 

83  *  INVERSE  USING  ADJOINT  AND  DETERMINANT 

-■ 

v  • 

84 

DET»2.0*AREAS<TRI> 

85 

MC  < 1 , 1 ) * ( X2*U3-X3*U2 ) /DET 

86 

MC  ( 1 , 2 ) *  <  U2— U3 ) /DET 

87 

MC(1,10>*(X3-X2)/DET 

88 

MC<2,1)=<X3*U1-X1*U3)/DET 

89 

MC  <  2 , 2 ) *  <  U3-U1 ) /DET 

90 

MC(2,3)-<X1-X3)/DET 

I 

91 

MC ( 3 , 1 ) » <  X 1 *U2-X2*U1 ) /DET 

92 

MC(3,2)*<U1-U2)/DET 

V* 

93 

MC ( 3 , 3 ) *  <  X2-X 1 ) /DET 

94 

95  *  ASSEMBLE  THE  NATURAL  COORDINATES  INTO  ARRAY  L<49,3> 

96 

DO  110  1*1,49 

L 

97 

DO  105  J*1 , 3 

98 

L  <  I ,  J )  *MC  ( J,  1 )  -h  MC(  J,2>#X<I>  +  MC<J,3)#U<I> 

99  105 

CONTINUE 

100  110 

CONTINUE 

101 

102  *  FIND  DELTA  U  AT  THE  SEVEN  LOCATIONS  WHERE  INTEGRATION 

1 

103  *  OVER  U  IS  NECESSARY 

104 

DO  120  1=1,7 

•  ’ 

105 

J»< I-1>#7+1 

.*« 

106 

DU  < I ) =U  <  J+6 ) -U  <  J ) 

107  120 

CONTINUE 

108 

r 

109  *  ASSEMBLE  DERIVATIVES  OF  NATURAL  COORDINATES  INTO  LX<3) 

110  *  AND  CALCULATE  INTERVAL  WIDTH  FOR  X  INTEGRATION 

■j 

111 

LX  < 1 )  =  <  U2-U3 > /DET  „  , 

G-3 


112  LX ( 2 )  =  <  U3-U1 ) /DET 

113  LX  <  3 )  =  <  U1-U2) /DET 

114  DX<TRI)=X<49)-X<1) 

115 

116  *  EVALUATE  M  AND  dM/dX  AT  THE  49  INTEGRATION  POINTS 

117  DO  150  1=1,49 

118  M<I,1)=L<I,1>**3 

119  MX<I,1)=3.0*<L<I,1)**2)#LX<1) 

120  M<I,2)=L<I,1)**2  *  L( I ,2) 

121  MX<I,2)=L<I,1)*L<I,2)*2.0*LX<1)  +  L<I,1)#*2  *  LX<2) 

122  M<I,3)=L<I,1)#*2  *  L< 1,3) 

123  MX<I,3)=L<I,1)*L(I,3)*2.0*LX<1)  -I-  L(I,1)*#2  *  LX(3) 

124  M<I,4)=L<I,2)*#3 

125  MX(I,4)=L<I,2)**2  *3»0#LX<2> 

126  M<I,5)=L<I,2)**2  #L<I,3) 


127  MX<I,5)=L<I,2)*2.0*L<I,3)*LX<2)  +  L<I,2>#*2  *LX<3) 

128  M<I,6)=L<I,2)#*2  *L<I,1) 

129  MX<I,6)=L(I,2)*2.0*LX<2)*L(I,1)  +  L<I,2)#*2  *LX<1> 

130  M<I,7)=L<I,3)**3 

131  MX<I,7)=3.0*L<I,3)**2  *  LX(3> 

132  M<I,8)»L<I,3)**2  #L(I,1) 


133  MX<  I , 8) =2  »0#LX  <  3) #L( I,3)*L(I,1)  +  L<I,3>**2  #LX(1) 

134  M<I,9)«L<I,3>**2  *L<I,2> 

135  MX < 1 ,9)“2»0*L< I ,3) *LX ( 3 ) *L( I ,2)  +  L<I,3)#*2  *LX<2) 

136  M<I,10)=L<I,1)*L<I,2)*L(I,3> 

137  MX  < I , 10) =LX <1)#L(I,2)#L(I,3)  +LX ( 2) *L< 1 , 1 ) *L< 1 ,3) 

138  MX  < I , 10) “MX (1,10)  +  LX(3)*L<I,1)*L(I,2) 

139  150  CONTINUE 

140  END 

141 
EOF  »  * 

EOT.. 


LI, lf 100 

1  *##*####**#####**#*#####*####*###*####**#*******#*##*#*#**# 

2 

3  *  THIS  SUBROUTINE  PERFORMS  WEDDLES  N=6  RULE  INTEGRATION  OYER 

4  *  A  TRIANGLE,  IN  THE  U  DIRECTION,  AT  THE  SEVEN 

5  *  X  COORDINATES,  FOR  USE  IN  THE  NUMERICAL  INTEGRATION  OVER 

6  *  SPACE  IN  SUBROUTINE  ISPACE 

7 

8  SUBROUTINE  ANING( DU ,M,GT, MX ,U,TRI , SIGMAT , SIGMAS, X ,CORDND ) 

9 

10  PARAMETER  (MNODE-151  ,  MNTRIA-46) 

11  DOUBLE  PRECISION  DU<7) ,M<49, 10) ,GT( 10, 10) , ILDF( 10, 10) 

12  DOUBLE  PRECISION  U<49> ,MX<49,10> , DFLUX ( 49 , 10 ) ,FLUX(49,10) 

13  DOUBLE  PRECISION  ML(MNTRIA, 10,10) ,MG<MNODE,MNOD£) 

14  DOUBLE  PRECISION  NLM< MNTRIA, 14 , 10, 10 ) ,LI < MNTRIA, 10,7) 

15  DOUBLE  PRECISION  NLI < MNTRIA, 7 , 10 ) , A3< MNODE*< MNODE-1 ) /2> 

16  DOUBLE  PRECISION  SIGMAT, SIGMAS, A, B 

17  DOUBLE  PRECISION  X  <  49  >  , CORDND ( MNODE , 2 ) 

18  INTEGER  TRI 

19  COMMON  MG, ML, NLM, LI, NLI, AS 

20 

21  #  CALCULATE  FLUX  AT  THE  FORTY  NINE  INTEGRATION  POINTS 

22  *  AND  CALCULATE  THE  DERIVATIVE  OF  FLUX  IN  THE  X  DIRECTION 

23  #  AT  THE  INTEGRATION  POINTS 

24  DO  100  1-1,49 

25  DO  50  J«l,10 

26  FLUX<I,J)»0.0 

27  DFLUX  < I , J ) -0  »  0 

28  DO  25  K-1,10 

29  FLUX ( I , J ) —FLUX ( I , J )  +  M< I ,K) *GT ( K, J) 

30  DFLUX  < I ,J) -DFLUX ( I , J)  +  MX ( I , K ) *GT ( K , J ) 

31  25  CONTINUE 

32  50  CONTINUE 

33  100  CONTINUE 

34 

35  *  CALCULATE  THE  INTEGRAL  OF  FLUX  OVER  U 

36  *  AT  THE  SEVEN  SPATIAL  INTEGRAL  POINTS 

37  #  PLACE  IN  ROWS,  THIS  IS  NON-LOCAL  INTEGRAL 

38  DO  200  1-1,7 

39  K=7*I-6 

40  DO  150  J»l,10 

41  NLI <  TRI , I , J) - ( DU< I ) /20  *0 ) * (FLUX ( K , J)+5  »0#FLUX (K+l , J ) 

42  C  +FLUX(K+2,J)+6.0*FLUX(K+3, J ) +FLUX < K+4, J) +5 ,0#FLUX <K+5 , J ) 

43  C  +FLUX<K+6, J) ) 

44  150  CONTINUE 

45  200  CONTINUE 


*  CALCULATE  INTEGRAL  U#DFLUX  -  ARANGE  INTO  COLUMNS,  ADD 

*  NLI  TO  OBTAIN  THE  LOCAL  INTEGRAL 

A— SIGMAS 

B».5*SIGMAS*SIGMAS  -  SIGMASfcSIGMAT 

DO  300  1-1,7 
K=7*I-6 
DO  250  J-1,10 


I 


56  ILDF  <IfJ)  =  <DU(I )/20.0>*<DFLUX(K* J)*U(K>+5»Q*DFLUX<K+1 , 

57  C  #U<K+l)+DFLUX<K+2f J)*U<K+2)+6.0#DFLUX<K+3, J>#U(K+3>+ 

58  C  DFLUX(K+4f J)*U(K+4)+5 .0#DFLUX(K+5, J)*U(K+5>+ 

59  C  DFLUX<K+6,.J)*U<K+6)  )#A 

60  250  CONTINUE 

61  300  CONTINUE 

62  00  400  1*1 f 7 

63  DO  350  J-lflO 

64  LI<TRIrJrI)*ILDF(I,  J)  +•  NLI  <  TRI  r  I ,  J>*B 

65  350  CONTINUE 

66  400  CONTINUE 

67 

68  END 

69 
OF  ♦ . 

EOT.  . 

UP 


» 


LI, 1,50 

1  ******************************************«*******************«X* 


3  *  INTEGRATE  OVER  SPACE  (X),  ACROSS  THE  LOCAL  TRIANGLE 

4 


5  SUBROUTINE  SPING<DX, TRI , TRIP ) 

6 

7  PARAMETER  <MN0DE=151  ,  MNTRIA=46) 

8  DOUBLE  PRECISION  NLK MNTRIA, 7, 10) , AS<MNODE*( MNODE-1 ) /2) 

9  DOUBLE  PRECISION  LKMNTRIA,  10,7) 

10  DOUBLE  PRECISION  UI 1 < 10, 10 ) ,UI2 ( 10, 10 ) ,UI3 ( 10, 10) 

11  DOUBLE  PRECISION  UI4< 10 , 10 ) ,DX ( MNTRIA ) ,NLM< MNTRIA, 14 , 10 , 10 ) 

12  DOUBLE  PRECISION  ML < MNTR I A, 10, 10) , MG < MNODE , MNODE ) 

13  DOUBLE  PRECISION  UI5 < 10,10 ) ,UI6< 10, 10 ) ,UI7< 10, 10) 

14  INTEGER  TRI, TRIP 

15  COMMON  MG,ML,NLM,LI , NLI , AS 

16 

17  *  TAKE  PRODUCT  OF  LI,  AND  INLF  -  RECALL  LI  IS  IN 

18  *  COLUMNS,  AND  INLF  IN  ROWS 

19  DO  100  1=1,10 

20  DO  50  J=1 , 10 


21 

22 

23 

24 

25 

26 

27 

28  50 

29  100 

30 

31 

32  *  DO 

33 

34 

35 

36 

37  150 

38  200 

39 

40 

EOF.  . 
EOT.  . 


UI1(I,J)=LI< TRI ,1,1) *NLI <TRIP,1,J) 
UI2< I , J)=LI (TRI , 1 , 2) #NLI (TRIP,2,J) 
UI3<I,J)=LI< TRI , I , 3 ) *NLI (TRIP,3 , J) 
UI4< I ,  J ) =LI <  TRI , 1 , 4 ) XNLI < TRIP,4 , J) 
UI5< I , J ) =LI < TRI , I , 5 ) #NLI ( TRIP ,5, J ) 
UI6 ( I, J) =LI <  TRI , I ,6 ) #NLI ( TRIP,6, J) 
UI7< I , J) =LI (TRI , I ,7) #NLI <TRIP,7,J) 
CONTINUE 
CONTINUE 


WEDDLES  N«6  RULE  INTEGRATION 
DO  200  1=1,10 

DO  150  J-1,10 

NLM<TRI,TRIP,I, J>»< DX ( TRI ) /20 .0 ) * ( UI 1 ( I , J)+5.0*UI2< I ,. 
C  HJI3<I, J)+6.0*UI4< I,J)+UI5( I, J)+5.0*UI6< I, J)+UI7( I , J) ) 
CONTINUE 
CONTINUE 

END 


■7 


Appendix  H  -  Spherical  Harmonic  Angular  Fluxes  and  Data  File 


SDATA 


This  appendix  contains  the  contents  of  three  data  files, 
PNDATA5  ,  and  PNDATA9  called  to  compare  finite  element  angular 
fluxes  in  subroutine  OUTPUT,  and  SDATA,  a  data  file  called  by 
subroutine  GDATA,  which  contains  submatrices  of  the  cubic  three 
dimensional  interpolate,  as  well  as  integrals  of  the  20 
polynomials  used  in  the  three  dimensional  cubic  fit. 

The  angular  fluxes  are  those  computed  with  46  legendre 
polynomials.  They  are  formatted  differently  in  this  appendix 
than  in  the  manner  the  code  of  appendix  A  reads  them. 


XE 

Pn  BENCHMARK  DATA  WITH  46  LEGENDRE  POLYNOMIALS 
FOR  c= . 5 


X 

U=-l  .0 

-.75 

o 

in 

♦ 

i 

-.25 

o 

o 

0.00 

0.5168E-01 

0 . 7383E-01 

0. 1009E+-00 

0. 1025E+00 

0. 1213E+00 

0.25  • 

0 .4860E-01 

0  ►  6105E-01 

0 . 7545E-01 

0 . 3841E-01 

0. 1028E+00 

0.50 

0.4154E-01 

0.4822E-01 

0 . 5639E-0 1 

0 . 6856E-01 

0.8217E-01 

0.75 

0 . 3335E-01 

0 . 3746E-01 

0 . 4289E-01 

0.5246E-01 

0 . 6354E-01 

1.00 

0.262SE-01 

0 . 2905E-01 

0 . 3295E-0 1 

0 . 4038E-01 

0 . 493 1 E-01 

1.25 

0 . 2050E-01 

0 . 2251E-01 

0 . 2544E-01 

0 . 3115E-01 

0 . 3938E-01 

1.50 

0. 1587E-01 

0  ►  1742E-01 

0 . 1972E-01 

0 . 2405E-01 

0 . 2986E-01 

1.75 

0 . 1223E-01 

0.1347E-01 

0 . 1532E-01 

0 . 1859E-01 

0.2320E-01 

2.00 

0 . 9401E-02 

0 . 1043E-01 

0 . 1 193E-01 

0 . 1439E-01 

0. 1803E-01 

2.25 

0.7221E-02 

0 . 8074E-02 

0 . 9299E-02 

0. 1116E-01 

0 . 1401E-01 

2.50 

0 . 5548E-02 

0 . 6259E-02 

0 . 7256E-02 

0 . 8665E-02 

0 . 1089E-01 

2.75 

0 . 4267E-02 

0 . 4857E-02 

0 . 5666E-02 

0 . 6737E-02 

0.8471E-02 

3.00 

0 . 3286E-02 

0 . 3773E-02 

0 . 4427E-02 

0 . 5244E-02 

0 . 6592E-02 

3 . 25 

0 . 2536E-02 

0 . 2934E-02 

0 . 3460E-02 

0 . 4086E-02 

0 . 5132E-02 

3.50 

0. 1960E-02 

0 . 2284E-02 

0 . 2705E-02 

0 . 3187E-02 

0 . 3998E-02 

3.75 

0 . 1519E-02 

0 . 1779E-02 

0 . 21 14E-02 

0.2488E-02 

0.3117E-02 

4.00 

0 . 1 179E-02 

0. 1388E-02 

0. 1653E-02 

0. 1943E-02 

0.2431E-02 

4.25 

0 . 9172E-03 

0 . 1084E-02 

0.1293E-02 

0 . 1519E-02 

0 . 1897E-02 

4.50 

0 . 7150E-03 

0.8466E-03 

0 . 101 IE-02 

0. 1188E-02 

0 . 1481E-02 

4.75 

0 . 5585E-03 

0 . 6620E-03 

0 . 7910E-03 

0 . 9292E-03 

0.1 157E-02 

5.00 

0 . 4370E-03 

0 . 5181E-03 

0.6188E-03 

0 . 7274E-03 

0.9044E-03 

X 

l)=  .25 

.50 

.75 

1.0 

0.00 

0 . 2755E+00 

0.5263E+00 

0 . 7671E+00 

0. 1014E+01 

0.25 

0 . 1821E+00 

0 . 3626E+00 

0 . 5886E+00 

0 . 8265E+00 

0.50 

0. 1300E+00 

0 . 2580E+00 

0 . 4498E+00 

0 . 6647E+00 

0.75 

0.9471E-01 

0 . 1863E+00 

0.3433E+00 

0 . 5328E+00 

1.00 

0 . 7029E-01 

0 . 1359E+00 

0.2619E+00 

0.4266E+00 

1.25 

0.5301E-01 

0. 1001E+00 

0 . 2000E+00 

0.3413E+00 

1.50 

0 . 4042E-01 

0 . 7435E-01 

0. 1529E+00 

0 . 2730E+00 

1.75 

0.3104E-01 

0.5561E-01 

0.11 69E+00 

0 . 2182E+00 

2.00 

0 . 2396E-01 

0 . 4185E-01 

0 . 8950E-0 1 

0. 1744E+00 

2.25 

0. 1856E-01 

0 . 3166E-01 

0 . 6861E-0 1 

0. 1393E+00 

2.50 

0. 1441E-01 

0 . 2406E-01 

0.5265E-01 

0. 11 12E+00 

2.75 

0 . 1 122E-01 

0 . 1835E-01 

0 .4045E-01 

0.8878E-01 

3.00 

0 . 8743E-02 

0. 1405E-01 

0.311 1 E-0 1 

0.7084E-01 

3.25 

0 . 6821E-02 

0 . 1079E-01 

0 . 2395E-0 1 

0.5649E-0. 

3.50 

0 . 5326E-02 

0.8310E-02 

0. 1845E-01 

0 . 4503E-0 1 

3.75 

0 . 4 160E-02 

0 . 6415E-02 

0. 1423E-01 

0 . 3588E-01 

4.00 

0 . 3252E-02 

0 . 4962E-02 

0. 1099E-01 

0 . 2857E-01 

4.25 

0 . 2542E-02 

0 .3846E-02 

0 . 8489E-02 

0 . 2275E-01 

4.50 

0. 1988E-02 

0 . 2985E-02 

0 . 6564E-02 

0 . 131UE-01 

4.75 

0. 1555E-02 

0.2321E-02 

0 . 5080E-02 

0 . 1440E-01 

5.00 

0. 1217E-02 

0. 1306E-02 

0 . 3934E-02 

0.  U44E-01 
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MICROCOPY  RESOLUTION  TEST  CHARI 

A 1  it  iNA*  I'ntuAU  t>5  *'•!  AND  ARPS  A 
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Pn  BENCHMARK  DATA  WITH  46  LEGENDRE  POLYNOMIALS 
FOR  c  =  »9 


s 


i 


X  U=-l  .0 


-.75 


-.50 


0.00 

0 . 1465E+00 

0 . 1566E+00 

0.25 

0. 1192E+00 

0 ♦ 1322E+00 

0.50 

0 . 1039E+00 

0.1153E+00 

0.75 

0 . 9149E— 01 

0 . 1009E+00 

1.00 

0.8035E-01 

0.8826E-01 

1.25 

0.7049E-01 

0.7725E-01 

1.50 

0.6181E-01 

0 .6763E-01 

1.75 

0.5417E-01 

0.5923E-01 

2.00 

0.4746E-01 

0.5188E-01 

2.25 

0 . 415SE-01 

0.4544E-01 

2.50 

0.3642E-01 

0.3981E-01 

2.75 

0.3191E-01 

0.3489E-01 

3.00 

0 . 2796E-01 

0.3058E-01 

3.25 

0.2450E-01 

0.2680E-01 

3.50 

0 . 2147E-01 

0.2349E-01 

3.75 

0 . 1882E-01 

0.2059E-01 

4,00 

0 . 1650E-01 

0.1805E-01 

4,25 

0.1446E-01 

0.1582E-01 

4.50 

0 . 1268E-01 

0 . 1387E-01 

4.75 

0.11 1 IE-01 

0.1216E-01 

5,00 

0 . 9745E-02 

0.1066E-01 

0 . 1708E+00 
0. 1484E+00 
0 ► 1289E+00 
0 » 1 121E+00 
0.9774E-01 
0.8538E-01 
0.7467E-01 
0.6535E-01 
0 .5723E-01 
0.5013E-01 
0.4393E-01 
0.3850E-01 
0.3374E-01 
0 .2958E-01 
0  *  2593E-01 
0  .2273E-01 
0.1993E-01 
0*  1747E-01 
0»  1532E-01 
0 » 1343E-01 
0.1178E-01 


0 . 1956E+00 

0 . 1636E+00 

0.1433E+00 

0 . 1252E+00 

0 . 1093E+00 

0.9549E-01 

0.8352E-01 

0.7309E-01 

0 » 6399E-01 

0.5603E-01 

0.4908E-01 

0 . 4301E-01 

0.3768E-01 

0.3303E-01 

0.2895E-01 

0 » 2538E-01 

0.2225E-01 

0.1950E-01 

0.1710E-01 

0 . 1499E-01 

0 . 1314E-01 


0.0 


0 . 2500E+00 

0 . 1915E+00 

0.1634E+00 

0 . 1421E+00 

0.1238E+00 

0 . 1081E+00 

0.9453E-01 

0.8273E-01 

0.7243E-01 

0.6343E-01 

0.5557E-01 

0.4868E-01 

0.4266E-01 

0.3738E-01 

0.3277E-01 

0 . 2872E-01 

0 » 2518E-01 

0.2207E-01 

0 . 1935E-01 

0.1696E-01 

0 » 1487E-01 


e 


i 


X  U= . 25 


,30  .75  1.0 


0,00 

0 . 3044E+00 

0.3292E+00 

0.25 

0.2396E+00 

0.2817E+00 

0,50 

0.1983E+00 

0.2387E+00 

0.75 

0 . 1690E+00 

0.2041E+00 

1.00 

0.1455E+00 

0 . 1757E+00 

1 ,25 

0.1261E+00 

0 . 1519E+00 

1,50 

0 . 1098E+00 

0 . 1318E+00 

1.75 

0 . 9575E-01 

0 . 1 146E+00 

2,00 

0.8366E-01 

0.9986E-01 

»  Am 

0.7318E-01 

0.8713E-01 

2.50 

0.6405E-01 

0 . 761 IE-01 

2.75 

0.5609E-01 

0.6654E-01 

3.00 

0.4914E-01 

0-5821E-01 

3.25 

0 . 4305E-01 

0.5094E-01 

3.50 

0 . 3773E-01 

0 . 4460E-01 

3.75 

0.3307E-01 

0.3907E-01 

4.00 

0.2899E-01 

0 . 3422E-01 

4.25 

0.2541E-01 

0 . 2999E-0 1 

4,50 

0 . 2228E-01 

0.2628E-01 

4.75 

0.1953E-01 

0 . 2303E-01 

5.00 

0 . 1712E-0 1 

0 . 2019E-01 

0.3434E+00 

0.3049E+00 

0.2690E+00 

0.2362E+00 

0.2068E+00 

0  » 1809E+00 

0 » 1582E+00 

0.1384E+00 

0 » 1210E+00 

0 , 1059E+00 

0.9261E-01 

0 . 8105E-01 

0.7093E-01 

0.6212E-01 

0.5440E-01 

0 . 4765E-01 

0.4174E-01 

0.3657E-01 

0.3204E-01 

0 . 2808E-01 

0.2461E-01 


0.3534E+00 

0.3228E+00 

0.2926E+00 

0.2623E+00 

0 . 2336E+00 

0 . 2074E+00 

0 . 1838E+00 

0 . 1626E+00 

0 . 1437E+00 

0 . 1268E+00 

0 . 1 1 1SE+00 

0 . 9856E-01 

0.8681E-01 

0.7642E-01 

0.6723E-01 

0.5913E-01 

0.5199E-01 

0.4569E-01 

0.4015E-01 

0.3527E-01 

0.3097E-01 
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